{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports & function definitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os \n",
    "os.chdir('/lustre/scratch/kiviaho/prostate_spatial/')\n",
    "\n",
    "import scanpy as sc\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import anndata as ad\n",
    "import squidpy as sq\n",
    "\n",
    "\n",
    "from scipy.stats import fisher_exact\n",
    "import matplotlib.pyplot as plt\n",
    "from scripts.utils import load_from_pickle, get_sample_ids_reorder, get_sample_crop_coords, get_sample_id_mask, save_to_pickle\n",
    "from statsmodels.stats.multitest import multipletests\n",
    "\n",
    "from itertools import combinations\n",
    "from scipy.stats import ranksums,kruskal\n",
    "\n",
    "import seaborn as sns\n",
    "sns.set_theme()\n",
    "\n",
    "sc.set_figure_params(figsize=(6,6))\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "samples = get_sample_ids_reorder()\n",
    "sample_crop_coord = get_sample_crop_coords()\n",
    "sample_id_masks = get_sample_id_mask()\n",
    "\n",
    "sample_categories_dict = {\n",
    "    'BPH':get_sample_ids_reorder('BPH'),\n",
    "    'TRNA':get_sample_ids_reorder('untreated'),\n",
    "    'NEADT':get_sample_ids_reorder(['goserelin','bicalutamide','degarelix','degarelix_apalutamide']),\n",
    "    'CRPC': get_sample_ids_reorder('CRPC')\n",
    "}\n",
    "\n",
    "region_colors_dict = {\n",
    " 'Tumor': '#fc8d62',\n",
    " 'Luminal epithelium': '#8da0cb',\n",
    " 'Basal epithelium': '#66c2a5',\n",
    " 'Club epithelium': '#ffd92f',\n",
    " 'Immune': '#a6d854',\n",
    " 'Endothelium': '#e78ac3',\n",
    " 'Fibroblast': '#e5c494',\n",
    " 'Muscle': '#b3b3b3'\n",
    " }\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Downloading ligand-receptor interactions results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_ligrec_results(results_dict,source_region, target_region, pval = 0.05,reverse=False):\n",
    "    # Function only reports those means that have a below-threshold p-value\n",
    "\n",
    "    if reverse:\n",
    "        source = source_region\n",
    "        source_region = target_region\n",
    "        target_region = source\n",
    "\n",
    "    for i,s in enumerate(list(results_dict.keys())):\n",
    "        means = results_dict[s]['means'][source_region][target_region]\n",
    "        pvals = results_dict[s]['pvalues'][source_region][target_region]\n",
    "        means = pd.DataFrame(means[pvals[pvals<pval].index])\n",
    "        means = means.rename(columns={target_region:s})\n",
    "        meta = results_dict[s]['metadata'].loc[means.index]\n",
    "\n",
    "        if i==0:\n",
    "            means_all = means\n",
    "            metas_all = meta\n",
    "        else:\n",
    "            means_all = pd.merge(means_all,means,how='outer',left_index=True,right_index=True)\n",
    "            metas_all = pd.concat([metas_all,meta],axis=0)\n",
    "\n",
    "    # Drop duplicate indices\n",
    "    metas_all.drop_duplicates(inplace=True)\n",
    "\n",
    "    return pd.DataFrame(means_all), metas_all\n",
    "\n",
    "def summarize_ligrec_results(ligrec_results,ligrec_meta,score_filter=10):\n",
    "    ## Summarize the ligand-receptor pair results\n",
    "    df_sum = ligrec_results.copy()\n",
    "    n_samples = df_sum.shape[1]\n",
    "    active_in = (~df_sum.isna()).sum(axis=1)\n",
    "    mean_activity = df_sum.mean(axis=1)\n",
    "    summarized_ligrec_results = pd.DataFrame({'mean':mean_activity,'active_in':active_in})\n",
    "\n",
    "    # Included n_references filtering\n",
    "    valid_meta = ligrec_meta[ligrec_meta['n_references']>= 3]\n",
    "    \n",
    "    summarized_ligrec_df = summarized_ligrec_results.loc[valid_meta.index.drop_duplicates()]\n",
    "\n",
    "    # Add columns for ligands and receptors\n",
    "    summarized_ligrec_df['ligand'] = summarized_ligrec_df.index.get_level_values(0)\n",
    "    summarized_ligrec_df['receptor'] = summarized_ligrec_df.index.get_level_values(1)\n",
    "\n",
    "    summarized_ligrec_df = summarized_ligrec_df.sort_values(summarized_ligrec_df.columns[1],ascending=False)\n",
    "\n",
    "    return summarized_ligrec_df\n",
    "\n",
    "def check_top_markers(r,p_thresh=0.05):\n",
    "    # These results rank the deg significance and specificity (low p-value = good marker)\n",
    "    deg_fishers_res_dict = load_from_pickle('./data/region_gene_markers_fishers_exact_test_240229.pkl')\n",
    "\n",
    "    fishers_test_df = deg_fishers_res_dict[r].rename(columns={'interaction':'gene','active_in':'deg in'})\n",
    "    fishers_test_df = fishers_test_df.sort_values('adj_pval').reset_index(drop=True)\n",
    "    fishers_test_df = fishers_test_df[fishers_test_df['adj_pval'] < p_thresh]\n",
    "\n",
    "    return fishers_test_df\n",
    "\n",
    "# Run the function straight off\n",
    "club_markers_df = check_top_markers('Club epithelium')\n",
    "marker_genes = club_markers_df['gene'].tolist()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ligand-receptor interaction region-specific prevalences (Supplementary Figure S8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Club region as source with all other regions as target"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "interact_with = ['Tumor','Luminal epithelium','Basal epithelium','Immune','Fibroblast','Muscle'] # Endothelium excluded due to n\n",
    "\n",
    "# Club as source in all samples\n",
    "source = 'Club epithelium'\n",
    "\n",
    "results_dict = {}\n",
    "n_interfaces = {}\n",
    "for target in interact_with:\n",
    "    dat_dict = load_from_pickle('./data/region_ligrec_analysis/{}_{}_slides_with_ligrec.pkl'.format(source,target))\n",
    "    \n",
    "    ligrec_results_out,ligrec_meta_out = get_ligrec_results(dat_dict,source,target,reverse=False)\n",
    "\n",
    "    summarized_df = summarize_ligrec_results(ligrec_results_out,ligrec_meta_out)\n",
    "    summarized_df['interaction'] = target\n",
    "    summarized_df['active_in'] = summarized_df['active_in']/len(dat_dict) # Define whether to normalize by the number samples with region interfaces\n",
    "\n",
    "    results_dict[target] = summarized_df\n",
    "    n_interfaces[target] = ligrec_results_out.shape[1]\n",
    "    \n",
    "club_as_source_ligrec = results_dict.copy()\n",
    "total_n_interfaces = sum(n_interfaces.values())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Format the data for plotting\n",
    "\n",
    "# Only include genes that are marker genes in the Club region\n",
    "df = pd.concat(club_as_source_ligrec)\n",
    "df = df[df['ligand'].isin(marker_genes)]\n",
    "\n",
    "ligands_to_include = df['ligand'].unique()\n",
    "\n",
    "df = df[df['ligand'].isin(ligands_to_include)]\n",
    "receptors = sorted(df['receptor'].unique())\n",
    "\n",
    "df['interaction'] = df['interaction'].astype('category').cat.set_categories(interact_with)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASAAAAEZCAYAAAA39vjlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAAsTAAALEwEAmpwYAABXIElEQVR4nO29d3wc9Z3//5yZrVpp1YstuTe5d5tqcCgG2wRD4tjGcFxCCBzgkHA54I4QSmJCwCFxyOUgPT9qIJhw9wVMCRCqjbGNu+WuLlld2r478/n9sdZKsopl491Zo8/z8fDDO/29o53XfMq7KEIIgUQikZiAarYBEolk4CIFSCKRmIYUIIlEYhpSgCQSiWlIAZJIJKYhBUgikZhGUgtQRUUF3/jGN7qs+/73v8/GjRv7PO7GG28kFArFza5zzz0XgB07drBmzZq4XQfgG9/4BhUVFSd93H/9139x5ZVXsnjxYu677z4MwwCgsbGR6667jksvvZTbbruNYDB4uk0+Zd5++21KS0vjeg3DMFi9ejULFizgsssu48MPPwTg3nvvZcGCBSxatIif//znJzzP2rVr2bJlS5/7rFu3jsbGxtNid2/8+c9/JhKJAFBbW8udd94Z1+sBbNu2jSuvvDL2b8KECezZsweI/r6+853vcPnll7No0aITf3+RxJSXl4ulS5d2Wfe9731PbNiwwSSLopxzzjkJu9bSpUtFeXn5SR/X1tYmhBDCMAyxatUq8fbbbwshhHj44YfFs88+G/v81FNP9fuckUjkpO04Ge666y7xz3/+s9/7G4YhdF0/qWs899xz4sc//rEQIvp9mpqahBBCfPjhh8IwDBEOh8W1114rNm7ceFLn7Ylrr71WHDhwoN/7n8r9nT9/vggEAid93OmioqJCzJ8/P7b8/e9/X7zxxhtCCCFaW1tFMBjs83hL3OUyjpx77rlcfvnlfPLJJwwdOpRf//rXaJrGV77yFV5//XXsdju//OUvef311yksLEQIwc0338zcuXM599xz+eijjwB4/PHHycnJYcWKFWzfvp2HH34Yv9/PkCFD+NnPfobT6ezx+hs3buT555/nF7/4BfX19Xz/+9+nqamJCy64gNdff5133nmHsrIy7rrrLgKBADabjYceeohRo0axbt06PvjgAxobG6msrOTWW2/lqquuQtd17rvvPjZv3syYMWNiLbna2lpuv/12/H4/hmHwy1/+klGjRvV6b1JTUwHQdZ1QKISiKAC8++67/O1vfwNgyZIlPPLII1x77bW9nufxxx+nqqqKQ4cOMX36dJYvX84DDzxAS0sLmZmZ/OxnPyMnJ4ctW7bw0EMPEQqFyMvL4/e//z319fXce++91NbW4nA4WL16NSNGjOC6665j/PjxbNiwAavVymOPPUZzczPvvPMOn332GS6Xi2eeeYbNmzezZs0ahBAsWLCAVatWAXDOOedwySWXsHnzZp588kl+9KMfcfToUQB++MMfMnfu3F6/zwsvvMBvf/tbADRNIyMjI/ZbArBYLIwbN47a2tpezwFw9913s3DhQubNm8dXvvIVlixZwttvv01aWhpPPvkkn3zyCTt37uS2224jMzOTZ599lvfee4/f/OY3BINBpk6dyv33309VVRW33norw4cPZ+/evbzxxhvcdNNN1NXVEQ6HWbVqFZdeemnM9qeeegpFUbj44ovJzs7m6NGjLF26lJEjR/KDH/yAO+64gxdeeAG/388999zD/v37cblc/PSnP2XEiBE8/vjj1NbWcvDgQerq6rj//vs577zz2LBhA6tXr0ZVVZxOJ88//3yf37+d9evXs2DBAgBaW1vZv38/jz32GABpaWknPkEcxfELc6IW0NixY8Wnn34qhBDipptuEh9++KEQouOtsG3bNrF06VIRCoVETU2NmDZtWuzYzq2YX/3qV+LZZ58VwWBQrFy5UjQ3NwshhPjd734nfv/733ezq/3YDRs2iO9973tCCCHuu+++WGviqaeeir0VfD5f7C2wbds2ccsttwghhHjppZfE4sWLhc/nE0ePHo3t/+qrr4rbbrtNGIYhdu3aJcaNGyfKy8vFH/7wB7F27VohhBChUEj4/X4hhBDf/va3RU1NTY/374477hBz5swRd9xxR6ylcN5558W2Nzc3i0WLFvX+Bzh2b6655hoRCoWEEEJ885vfFBUVFUIIIV577TXxwAMPiGAwKC6++OLY2769VXHHHXeInTt3xr77d77zHSFEtGXw0EMPCSGEeOutt2LrO7eA/H6/mD9/vqiqqhKhUEgsW7ZMbNq0SQgR/bu///77Qggh1q9fL+68804hhBC6rsdafv/1X/8ltm/f3u37nHfeeeKxxx4TV111lbjjjjtES0tLl+0ej0dccskloqqqqs/70tnW+fPni5dfflkIIcSDDz4o/vrXv8a+Z/s9aWhoEP/6r/8aa63cf//94vXXXxfl5eViwoQJYt++fbFzt9+/trY2sXDhQmEYhtizZ4/46le/Gvt+7ft0bgF1fl5++9vfxlp67733nrj++uuFENG/5ze/+U0RiUTEjh07xLJly4QQ0efnk08+EUJEWy5CCFFTUyO+/e1v93kfvv71r4tt27YJIYTYtWuXWLlypbj99tvFlVdeKR577LE+jxUiyVtA7W/t3ta73W5mz54NwPjx46msrOyy35YtW7jkkkuwWq3k5+cza9asPq93+PBhSkpK+Jd/+RcAwuEwZ599dr9s3bp1K9/97ncBWLhwIX/84x8BCIVCPPjgg5SUlKCqapcxl7PPPhun04nT6cQwDMLhMFu3buXyyy9HURQmTJjAyJEjAZg8eTJ33303qqpy2WWXMXr0aAB+97vf9WrTz3/+c0KhEPfccw+ffPJJ7C1/slx88cVYrVY8Hg9btmzhlltuAaLjKUOGDOHQoUMMHTo01iJrb1Vs2LCBAwcOxM6jaVrs86JFi2LnfuCBB7pd8/Dhw4waNYpBgwYB0Xu6ZcsWZs2ahcvl4vzzzwdg7NixPPzwwzz66KMsWLCAKVOmALB69eoev4vH46GoqIh169bxxBNP8N///d/853/+JwBCCO655x6WLVsWu25/ueiii4Do77C8vLzb9s8//5ySkpLYmGYgEGDw4MFMmjSJkSNHMmbMmNi+f/7zn3nnnXcAqKyspK6ujk8//ZRFixbFWrbt97g3tm7dys033wzABRdcwA9/+MPYtgsuuABN05gwYULsmZk+fTqPPvooV111FZdffjkA+fn5ff6+KisraWxsjN1zXdf5/PPPefnllxk+fDi33nor//jHP2L3pieSWoDS09NpaWnpsq696Q9gs9li61VVRdf1bufoLGK9fW7v5gghmDx5ckw8TgbRS0jdX/7yF0aMGMGaNWtoamri61//emzb8fa3DxT3ZOfs2bN55plneOedd1i1ahU/+tGP+iWONpuNSy65hLfffptzzz2XlJQUPB4Pqamp1NTUkJeXd8JzOByO2HfMz8/nlVde6bJ97969PR6nKAovv/wyqtr3XEdvL5re6NwlHjFiBC+99BLvvfceDz74ICtWrOBrX/tar8fm5eXFHohLLrmEBx98MLbtV7/6FTabjRtuuOGk7IGOv2Vvv0MhBBdddBE//vGPu6yvqKjo8n02bNjAjh07+Nvf/obNZmPx4sWnfUKls63tv7mbbrqJefPm8c4777B06VJeeuml2HPWG2+88Uas+wXRezt8+PCYmF544YXs3bu3TwFK6lmw1NRUnE4n27ZtA6J/rPLycoYNG9av42fMmMHbb79NOBymtraWzz77LLbN6XRSU1NDKBSKjQWNHDmSiooKSkpKAPD5fP2elZk+fTqvv/46QOx/iL5xc3JyYg/jyZxn7969HDp0CIi+bXJzc7nmmmtYtGhRzMaeEEJQVlYGRFsp7777bqwldeGFF/J///d/APz9739n/vz5ALz11lsnnP1JS0vD7XbHZo7C4TAHDx5k5MiRlJeXc/DgQQCam5sBmDlzJi+++GLMjs42t3/Hd999lwkTJgDgcrnwer1AVFgOHjxIbW0tkUiE9evXM2PGjG421dbW4nK5uPrqq1m5cmWvYtjOhRdeyKZNm4DoGF77fXnllVf49NNP+clPftJl/6effpqnn366z3P2RufvM23aND755BNqamoAaGpqin3ujMfjIT09HZvNxvbt22P39KyzzuLVV1/F4/EAHfe48zU6M336dF577TUA3n///ViLuTfKy8sZP348t956K4WFhT3adjzr16+PtZYg2mJKS0ujuroaIQSbNm2K3d/eSOoWEMDDDz/Mgw8+iNfrxWq1snr16i4th76YMmUKZ511FosXL6aoqIhp06bFtq1atYprr72WgoICRowYAUTfDGvWrOG+++7D5/MhhOA//uM/+iV4t912G9///vd5/vnnOfvss3G5XAAsX76c22+/neeee67PN0E7CxYs4KOPPuLyyy9nzJgxjBs3DoBPP/2U3//+91itVjIyMmIDfTfeeCM/+clPyM/Pj51DCMGdd96J1+tFCMHMmTNZvnw5EH3Tffe73+WPf/wjY8aM4fbbbweiP8D25n1frFmzhvvvv59HHnkEXde54YYbGDVqFD/72c+48847CYfDDBo0iCeffJJ7772X++67j2eeeYZIJMKSJUti30fXdb761a/GBqEh2s269957eeKJJ3jmmWe47777uPHGGzEMgwULFvTYhd63bx+PPPJIbPD0Zz/7GQD33HMPy5cvZ/LkyV32v+mmm7jjjjv49a9/TW5uLo888ggADz30EGlpaSxduhSAb37zmyxZsoTDhw8zffr0E96Xnrjqqqu48847ycrK4tlnn+VHP/oRt9xyC+FwGKvVyo9//GPS09O7HHP++efz7LPPsmjRIsaNG0dxcTEQ7Wpee+21LF++HFVVueSSS1i1ahVLly7lmmuuobi4mB/84Aex81x77bXcc889XHHFFbFB6L7405/+xMaNG1FVlRkzZlBcXExtbS0//OEPe+yGVVVV0djY2O3+3n333fzbv/0buq4zffr0Li2kHjnhKNGXiHhO4QcCgdg06t///ndx1113xeU68eLOO+8UDQ0NCbnWyU5Pm8nNN998wqlkyamT9C2gM4WKigr+/d//HV3XycjIiL2JzxTONHsTxf/8z/+YbcKXGkUImZBMIpGYQ1IPQkskki83UoAkEolpnHFjQIZhxGbETtZ/RCKRnH6EEITDYVwu1wl9vo7njBMgr9fLvn37zDZDIpEcx9ixY/sX/9WJM06ArFYrEP2y/fUHkkgk8SMUCrFv377Ys3kynHEC1N7tstls2O12k62RSCTtnMqQiByElkgkpiEFSCKRmIYUIIlEYhpSgCQSiWnEVYDC4TDLly9n1qxZrF+/vtv2d999l2XLlrF8+XK2b98eT1MkEkkSEtdZMIvFwq9+9Sv++te/dtum6zpr167l6aefxuv18r3vfY/nnnsunuZIJJIkI64CpChKrxn3jhw5wvDhw0lNTSU1NZVIJEIwGJRT6xLJAMI0P6CWlhbcbnds2e1209zc3CWxVl/s3LkzXqZJJJIEYZoApaen09bWFltua2s7YaLtzkyaNEm2liSSJCAYDJ5yg8A0ARo2bBhHjhzB5/Ph9XrRNE0KikQywIi7AN1+++3s3LmTlJQUtm/fTkZGBhdffDEjR47ktttu45vf/CaKosRKo0gkkoHDGZcRsb25J7tgEkly8EWeSemIKJFITEMKkEQiMQ0pQBKJxDSkAEkkEtOQAiSRSEzjjMuIKJGcqfjrD+A7WoLVlUXakDkoqma2SaYjBUgiSQBhbwNNJW8CgmDTERTVStqQ7rXuBxqyCyaRJAA92AZ0uNxFAi3mGZNESAGSSBKALX0wFlc2AIpqwZU/wWSLkgPZBZNIEoCq2cid/HXC3jo0hxvN5jLbpKRACpBEkiAUzYLNPchsM5IK2QWTSCSmIQVIIpGYhhQgiURiGlKAJBKJaUgBkkgkpiEFSCKRmIachpdIBiiNjY2xgqCTJk0iJycn4TbIFpBEMkDZunUrHo8Hj8fD1q1bTbFBCpBEMkCJRCI9fk4kUoAkkgHKpEmTUFUVRVGYNGmSKTbIMSCJZIBSWFjIoEHR0BBVNactIgVIIhnAmCU8seubenWJRDKgkQIkkUhMQwqQRCIxDSlAEonENKQASSQS05ACJJFITEMKkEQSB4QRIRJoRQijy3o95EMP+02yKvmQfkASyWkm4m+hfufLGCEv1rQCciZeiaJZ8FRupfXIx4BC+sh5uAaZ432cTMgWkERymvFW78AIeQEIt9UQaDyMEILW0o3H9hC0lm0wz8AkQgqQRHKaUa2O45adKIrSZb1mTUm0WUmJ7IJJJKeZ1MLp6ME2wp46HDmjsWcUAZBVfDmtRz4i2gU731wjj2EYBrt27aKxsZH8/HyKi4sTev24C9ALL7zAunXrsFqtPPTQQwwZMiS27eWXX+bpp59GVVUWL17M9ddfH29zJJK4o6gaGaPnd1tvS8snZ/LVJljUO6WlpZSWlgLQ1tZGRkYGBQUFCbt+XLtgzc3NvPjiizz99NP8x3/8B2vWrOmy/cknn+Spp57ir3/9K3/9618JhULxNEcikRzH8c9cop/BuArQ9u3bmTNnDhaLhSlTpnD48OEu20eOHInP5yMQCOBwONA0LZ7mSCSS4xg2bBgpKdHxqPT0dAYPHpzQ68e1C9bS0kJ6enpsWQjRZfuCBQtYsmQJmqZxww03nJQA7dy587TZKZEMZNxuN6mpqaiqyrZt2xJ67bgKkNvtpqSkJLbcOfeIx+PhySefZP369dhsNr71rW9x8cUX91uBJ02ahN1uP+02SySSkyMYDJ5ygyCuXbCpU6eyadMmdF1n165dDBs2rOPCqorVaiUlJQWbzYbD4cDj8cTTHIlEkmTEtQWUkZHBkiVLWLlyJRaLhdWrV7Nu3TqKioqYM2cOV111FcuWLUNRFKZNm8bYsWPjaY5EIkkyFHH8wEyS097ck10wiSQ5+CLPpHRElEjihKGHaDnwHmFvPc68caQVzTTbpKRDCpBEchrYu3cvu3fv7rJOD7Ti9UbHNR3WT7C69qNottj2CRMmJNzzONmQsWASSZwQwiAQhkC4Y1nSFdkCSnK27qml5EgT+dkpzJtVhGZyGRVJzxQXF3drzUT8zfzthedA6FwyexhZ4xeiqNLZtjNSgJKYipo23v20HICqox5SU6zMnjTIZKsk/cXizMCWlg/CIGvCYhRFMdukpEO+TpMYjz/cddkX7mVPSVKjqFJ8ekEKUBIzsiidvKxonE6Kw8LUcbkmWySRnF5kFyyJsVk1ViwsprktSJrLhs0qxw8kp5/W1lb27t2LoiiMHz+e1NTUhF1bClCSo2kq2RlOs82QfInZtGkTfn80Ub7P5+OCCy5I2LVlF0wiGcAIIQgEArHlzp8TgRQgiWQAoygKo0ePji13/pwIZBdMcsaiGyEEAosqYwK/COPGjaOoqAhFUWLJyRKFFCDJGcmR1g/ZWPPfCGEwM/8GxmRcarZJZzQul8uU60oBkpxRtMdc1fv3YTARgHfYxEZxFKDHN7iMuUpe5BiQ5IxEUTp+ugoqPp8Pn89nokWSU0G2gCRnFO0xVw3+/Xxa+1uE0JmZ/y0+Wr8PgKuvTq6yN8lKY2MjBw4cwG63M378eGw224kPigNSgCRnJNnOMVw+/NFOa/aZZsuZRjgc5tNPPyUSiQAQiUSYOdOcXEWyCyaRDDDC4XBMfICYE6IZSAGSSAYYTqczVv1UURRGjBhhmi2yCyaRDDAURWHmzJm0tLRgs9kS7vvTGSlAEskARFEUMjIyzDZDdsEkEol5SAGSSCSmIbtgEskAxDAM9u/fT1tbG4WFhQwaZE6qX9kCkkgGIAcPHmT//v3U1NSwZcsW2traTLFDCpBEMgDxer2xz0II08JYpABJJAOQIUOGoB4r8ZSWlkZ2drYpdsgxIIlkAJKdnc38+fPx+Xykp6djsZgjBVKAJJIBitPpxOk0N9+47IJJJBLTkAIkkZxm9GAbhh6KLggj+k/SI7ILJpGcJoQQNO97C3/9fhTNhqtgIqG2GgDaKraQVjTDZAuhoaGBnTt3IoRgwoQJ5OXlcfjwYZqamigoKGDw4MEJtUcKkERymoh46/HX7wdA6CE8Vdti29rKNpJaOK1LJsdE0Z7GVghBU1MTQgh0XWfv3r243e4uU/Lp6elMmTIlYSlspQBJJKcJxWIHFEBEl1UroAOgWuymiE9nDMNAiKhtuh61q3NeoM7rE0XcBeiFF15g3bp1WK1WHnroIYYMGRLb1tDQwAMPPEBzczO5ubn8/Oc/j7c5EkncsDjcZIy5CG/VNjRHGq5BU1D2vwpAVvFC0+xqT2NrGAbvv/8+Ho+H6upqNE1j8eLFbNy4EV3XcTgcnHfeeTgcjoTZFlcBam5u5sUXX+S5555j9+7drFmzhrVr18a2P/zww/zgBz9g6NCh8TRDIkkYKXnjSMkbF1u2unIAsLkLzDIphqqqnH/++VRWVtLW1obdbicrK4sLL7yQtrY2MjIyEp4bOq5twu3btzNnzhwsFgtTpkzh8OHDsW26rnPo0CEef/xxrr32Wl577bV4miKRSABN0xg6dCh2e0cxR6fTSV5enimJ6ePaAmppaSE9PT223N7/hGj3q6SkhEcffZSCggKuueYazjnnnH4nSdq5c+fpNldyBuPxeADYvHmzyZYAQpAR3o/DaCLQqhBRUpLDrk4ky/2KqwC53W5KSkpiy+2xJxAdbR88eDAjR44EYOLEiZSVlfVbgCZNmtRFxc90DlU089bHpQghuPTc4YwsyjDbpKRGN8K8duT7eMK1WFUXTtflqIpGS96rNPj3UZg6h7MG3YqqaAm3zVe3j+Z9jQCoGDhtmFZ1ojdKS0uB02NXMBg85QbBSQlQWVkZlZWVXUbKzzvvvF73nzp1Kr/5zW9iU37Dhg2LbbPb7eTn51NfX09mZib79+9PuA9CMvHmR0fwBaIzEm98dIR/WzbNXIOSnO31z+IJ1wIQNrz4QpXYVBch3w4ASts+oCh1FkPd5yTcNgUl4dc8VYLBIAcOHEAIwejRoxM6AA0nIUAPPPAAW7duZdy4cV1aMn0JUEZGBkuWLGHlypVYLBZWr17NunXrKCoqYs6cOdx5553cfvvthMNhrrjiCnJycr7YtzmDEZ0/i153kxyjP7dI9Guv048jZxTOprEEGg6haCqaLdUUO/rD5s2baWyMttYaGxuZN29eQq/fbwHasGEDr776ahfx6Q8rVqxgxYoVseXOraDJkyfzzDPPnNT5vqxcevZw3vj4MELApecMO/EBA5ypOSsob9uAL1KPRXHitg1GVSykO8fTENhPYepshqSdZYptiqKSOfYSAKzV60yxob+0trbGPpuRlKzfAjRx4kQqKyu7+PFITh+jhmZwy9DpZptxxqCpNq4c9T+x5XU7ow/6xUMfNMukM5KioiKOHDkS+5xo+i1Ahw4dYvHixYwaNQqbzYYQAkVReP755+Np34DDMAQtbUFcKVZs1sQPoEriQ/2OdYTaatHsaeRM+TqaNbFjLb0xadIkCgoKEEKYMgTSbwF6/PHH42mHBAhHDF56s4SqOi9Oh4Wll44jJ9PcfC2SL44ebCPUGp321gMtNO17i5yJV5hsVQdmjr2eUIAikQgWi4Xc3NxE2DOgKa1qoaouGhjoD0T4fO9RLj5bjged6Qija3yViARMsiT5OKEA3XrrrTz55JNcdtllKErH9GJ7F+wf//hHXA0cSLic1j6XJR34I020hqrIsA2nLVxBRASxKHY8oVp8kQayHWPQ1OS4fxaHG0WLIPQwoOAe0fvM8UDjhAL05JNPAvDOO+/0ud/nn3/OtGnTTotRA5VBual8Ze5Q9hxqIDczhdmTzY8fSkaaAof5R/l9hA0/FsVBRARoCUzFqWXw6uHbMdDJcYzjK0PvQ1OSQIQUlfw53ybib8DiyEDVksCmJOG0xYI98MADp+tUA5ppxXmsWDiei88ehkWTCSt74nDr+4QNPwAR0dGdCegtGMfSX9QHSmgKHDHDvB5RVRWbK1eKz3GctlAMIb3nJAkizdZzFU9V6fg5a4qVFEtWokw6Ib66fQTq92N15ZI6ZJbpuYGShdMmQJ3HhySSeDI6/RLCuo+GwAGyHaNoCpZyQLOQYs1ldOZgvOFaRmcsIMVqTq2r4xF6mOZ9bwEQaDyCotlILZxmrlFJgmwBSc44FEVhQvaSLutqbVFHxOl5y0ywqG+E0TXrYCTQYpIlycdpawfOmGF+wm2JJBlRLQ4szkwAFM1KSt54ky1KHvotQKWlpdx88818/etfB2Dfvn387ne/i23/0Y9+dPqtk0i+DCgKOVO/Tvbkq8ibsRJbWp7ZFiUN/Rage+65h1tvvTWWxHrMmDG88sorcTNMIvkyoWo27O7BaDaX2aYkFf0eAwoEAkyePDm2rCgKmiZjlb7MhA2ddypL8OkhLhw0lkx7itkmdSFiBClpeo2I4ccQli6zYGYT9tThq9uHEfKh2pLrvgUCAXw+H4qioOu6qc9xv/9ieXl57NmzJzbb9cILL8jI+C8xe/fu5Z3NG2gNRf1sDqubGJ6Wjc/nAyAlpftDNWHChITVkwLYWPPflLV9AoA3OIcsx6iEXbsv9KCH+p0vI/QwkYBAI7kqo27YsAG/P+pH9fnnn5uarbHfXbAf//jH/P73v6euro7zzz+fTz75RDoffskJRMKxzyFDRxcGPp8vJkJm0+A/EPscEUFEkjzoEX/TsbCLKKLTfTSbSCQSywcN0co1ZtLvFlBaWpqs2zWAKC4uZq8jyOvluwAY7c7l61Mu5uWXXwbg6quvNtM8AIrS5lDSFK27ZVNdKPEt8tJvrKm5aPZU9GD0QVeTJPUGgMViIScnh+rqagAKCswN9+m3AC1evJiJEyeyYMECLrzwwoTnjpUkniXDpzLKnYM3EmJG9pCkczadnns9eSkTiRgBttiPmm1ODNXiIGfKUoJNpVirtqFoiS930xezZ8/m8OHDUX+qCRNMtaXfAvTGG2/w2Wef8cYbb/DYY49RXFzMZZddxsKF5lV8lMSfyVmFZpvQK4qiUJQ6G4AtJFfqU82WQkr+eBRtj9mmdEPTtFhFGbNfKv1usyqKwuzZs/nhD3/IU089hdVq5d///d/jaZtEckYT9jVghLzdPKElHfS7BVRTU8Obb77Jm2++ic/nY8GCBaxatSqetkmSDE84gIAzqOiMeYQ9ddTteInIsVJLIU8dtlSZ1O94+i1At9xyC5dddhkPPfTQl6KWu9AjiG3vQiiAMvVCFGea2SYlLUII/ljyMZ/WlTK61UORK8Nsk5KeYHM5dMqEGGwuSyoBCofDRCIRampqaG1txWazMXTo0JOuevNF6bcArVuXXH3sk2Xv3r3s3r27Y0VjDfha8ClW2HyQlEFdRTXRPi3JzBFPA5/WRStp6sKgMeg12SIoa/uET2ueiOUFavFPAeCFfS+RaR9OyPAyKv0SirMWmWKfNa2AaFsxGqRtS0ue5HK1tbWxcjyfffZZbH1bW1sXZ+NEcEIBuvvuu3n44YdZtmxZjylZz9iqGKHoD9eHBQxIMQxIsPqfKTg1a6dHCVSTc9lEjAAfV61FoHfbposQ9YF9AGyt+zP5KRPIdIxItInY0weTPfEKtNJ/oFjs2NOTZzC/qampx/XtBQoTyQkF6I477gDgsccei7sx8aS4uLhLi8Z451nE5//gFctosDm5+liQraQ7BSnprBg1m3erSki1BshxmBvPpItwj+LTEyHDPKdJe8YQNIfbtOv3Rm5uLoqidEuhk5+fn3BbTvgqy8uLRu7+/Oc/p7CwsMu/M9kxUZm/AmXRTZBZALkypOREXDB4DPfPWszglAw0k1tAdi2NMRmX97pdU6JTzEWpc8lzmp/6QhgRgs0VGJEQweYKIv5mU+3Jzs7G7Xbjcrk499xzmTx5MjNnzjRlyOGELSDDMIhEIhw6dIhwOBxTTY/Hw549yefj0F8URUEZNwd2VZhtiuQUmJX/LSZlf42IEUQIwRs7/4khInx11HexaSmEdC8OS7rZZmKE/UT8TTTsegVFsyH0EKCQOe5SnDmjTbPLYrFgsVjIzMwkMzPTPDtOtMNTTz3FX/7yF44ePcpll10WE6DU1FSuueaauBsokfRGZ4FRlWg0vMPi7rbNTPRQx4B9VHwABL7a3aYKULJwQgG6/vrruf7663nmmWdYuXJlImySJBE1vlbeqtyDQ7Pi0CzU+FtJt8pqrf1F0SydhKcDzZGReGOSkH5Pwzc1NdHa2orbHX3DtLS08NRTT3HbbbfFzbh4IJqPIko2oWTmoYydbbY5SY1uGPxixz9oPjZjCDAiFMATDlIf8JDjSDXRujMDiyMdHRVHziBsaYMJtVWh2dNIG5LY356u65SWRl0phg1Lnmq7/Ragt99+u4vYpKend1uX7Ah/G8ZzD4G/LerRe2HyJwdv8QSxWTWc9sQn2/JFQl3Epx1DCOr8ySNAQhjoIoSqWAjqXpoCB8l2jMaqJUMiMAXN4SZr3IJjy4n1s2nns88+o66uDiD2fzLQ71+1ruv4fL5YIiqPx4Ou928qNGmorwJ/W2xRlO8FBptnzwl446Mj7DpQj6YpXHHhKEYWZST0+mk2B5OzBrOjsQroCMGwqxoj0pKj5I0hdD6o/BkNAQMFlZcPvIjAQMXCZcMfJd1elHCb9JCPsLcBa8qxumTCINRaizUt17R6YA0NDT1+Npt+C9B1113HypUrWbhwIUIIXnvtNa6//vp42nb6ySkEVzp4oy0fZfhEONizU1ai6OahfQyPx0tLWxDLsfrmb67fQ15WSsI9tP9twjz2NteQotlwWKy8Xf06Ts2Kw5IcFT4b/Puo8m4FpiIwYknJDCJsq3uWeUV3JtQeT9U2Wg9/GFuO+B0YYT/1O/6GYrGTP/NaVEviU9nk5uZSW1sb+1xeXp5wG3qi33L8jW98g5/+9KeEQiFSUlJYs2YN55xzzgmPe+GFF1i+fDnXXXddj1+6ra2NuXPnsn79+pOz/BRQnKmoK+5BuWAZ6pWrUKfOj/s1TxV/wA+iI5OeppoTAqopKhMzBzPCncOglHRcFhtqEuUFslvcKL2Exzqtia+M2lq6ocuyEe7owopIEE/1jkSbBETLZk2ePJnJkycnVQmtkxoD+uUvf0llZSVFRUX89Kc/Zdy4cX3GiDU3N/Piiy/y3HPPsXv3btasWcPatWu77POHP/yBqVOnnvo3OEkUdzbKzEsTdr0TcbyHdjvr1q0jEIygpE/G5bTylblDSU1JrsRWyYDbVsicglt4V92MptjIso+iJVRBpn04M/O+mXB7VIsdI9R7+g3Nas64lKZpSTX43E6/BehXv/oVzz//PCtXruSVV15h+/bt/OUvf+nzmO3btzNnzhwsFgtTpkzh8OHDXbbX19dTXl6e8AC4MwWH3cLVC8335E12RqZfyOf2aBzTguHmpojJGr+Qxt2vYoR9KBY7lpQUdH8LqAr2jKGk5JubgTDZ6LcAWa1WUlOjsx6RSIQpU6awb9++Po9paWkhPb3DIez42JMnnniCG2+8kbfeeutkbAZg586dJ31MT7Qn6N68efNpOd/pot2uTZs+o6ENrBpkpJrf9Un2+5UUdmlT4FilG19gFyguqm0TwQdHtmwx17ZjJMv96rcA5ebm0trayvz587n55ptJT0+PxYn1htvtpqSkJLbcOddIeXk5ra2tFBcXn5IATZo0KZZW8otQWloKQjBdawGLDWXC2Siq+fXO2n02jjSmUV4Tnbk7d/og5k4xd9au3S4zS7n0RCLt+uCDD/o9lR0IRMsatdvXH3Jzczn//PNPybb+cuTIEXw+H6qqMnbsWFyuUw8wDgaDp9wg6LcAPfHEEwB873vfY+PGjXg8nhPepKlTp/Kb3/wGXdfZu3dvlz7onj17KCsr44YbbqCsrAyXy8WoUaMYM2bMKX2RL0R9BaL0tejnihKUy25IvA09oBsiJj4Aew83mi5AieRkHvT6+nrg5PJWneqDXldXR319PVlZXQe5NcVgYlYTLmvHGFCZRSAA1dZKhefEflOJSIlhGAbNzc0YhkFlZSXV1dVcdNFFp+WFfrKcknfb3Llz+7VfRkYGS5YsYeXKlVgsFlavXs26desoKiri0ksv5dJLo4PBjz/+OGPGjDFHfISAYEfKBlG6K/E29IKmKqS5bLR5o678eVnJ4FiXOHp70Nsx3OUIawC1eShOZzQ8JJJSibC3obYWooR7f6t/0Qc9KyurW0GGQM02AtVdBXNYdnu32cvEuZeh2ftOz/Haa699Ibv6QzAYxDA6aqgZhkFLS8sJezTxIO7utStWrGDFihWx5Z5G4k3NLa0oYHPCsZeWUpRcWRCXLhjHlt21OOwasyclT1a9RNHTgw6wt2UdJa0fA5AytIbLC35MufdDPm/6XwCsRYf5SsFqbFrPrY5EPOjdMX8MD8But6OqakyEVFWNhVglmuQppm0muUUow8aDxYoyeZ7Z1nQhI83OV+ae+Tm4Tze1gQ5/Gk+kCl+kocu6oNFMS7icXC1xs4j23AmEW8rQffXQJYekgmPQDDR7cuQdV1WVjIwMvF4vBQUFjBs3zrQ6f1KAABQVddaCE++XQISANm+Iv71ZQl52CrMmFpDiSA7v42Qgx15Mc+gQAClaLk5LJjn2Yqr90RzHNjUVtzWxaVAVzUrauCsSes1TRVEUUlNTmTVrlql2SAFKUhqa/QRCERqr2yirbuNIZSv/8tWJptnzbtU+djZW4Qz6sGsWfr3rPXIcqVw1fBp2LfE/ownpS0m1DCJotDDEdR6aYmNk2iXYNTeecA2FKXOxa8mXDlXSFSlASUoo3DXQt77JTziiY7Uk3kVgR2Mlzx+MtixGBDwowKHGjkHM5aMS/xZVFJVhqd27y4Up/ZsgiSdG2I+v7COMkAd73kTs2dHJlXBrJYGqzaBZSRlyDpojOZKmmcmAFSCj+hDio5ehTgVX8r0pHXYLvkBHLNiwQW5TxAegPtC1DI/oss2TWGPOAPxVm4i0RuMe/WUfYUktQLWm4D38Dhyrkuor+4C0sYsTalcoFKK6ujo23hOJRNi+fTvhcJisrCyGDh2KpiX2NzYgBcg4WoZ4bnV0wTIamgMY7zyD+pXkyfiYle7AYbdQODqP3MwUikcmPrCynRk5Q3irYg8NQS+aomI75qhpUVQuGGSC60SSIyKdMyCKaEZEiyMmPt33iT+6rvPRRx/h9UZfJh6Ph2AwSFlZGQDV1dXU1NRw9tlnJ9SuASlAHPy82ypxaBskkQABpDgszJ9j/gxYus3JvTMWUuVrZuvR91EVha9OnUe63UmW3dwSPcmIo2AKXm8tQg9hzRyF5sxCURTsBdMI1nwOioZjcGIj0r1eb0x8IFoZ9XgaGhrQdT2hraCBKUBDiuGTV7quG5RcCcJ1Q9DcFuRPL+/EZlUZPTSTuVMGmWaP02JllDuXbcdScYxw55hmC4AnXMOO5mdoDVWAIrApLoJGGyHDg11N59y8u0i1muM3ZXHl4Z60DKGHUDtFvzsHTceeOx5F0VC0xM5opqSkYLfbCQaDURstFkKhrq0wt9stu2CJQC0ai7H436JjQF4LpKSiLvy22WZ1obElQDAUIdIajSWqbfCRnmqjeGRyZCI0m0/r19IWqYotB+hILBcwGvno6MMsKPylCZZFUVQLitr98TIjGRlEBeecc86hoqICh8PB5s2bCYVC5OTkEIlEyMzMNCUSYUAKEIA6dhaMnQXHYofMSpXZG7pudFvX6k3suAHArsZqKr1NFKZmMNTVdRzqcGs9Vk2jyJX4ulI+ve9QipBhzuC4EfYRqNqMHmhFUTU0Vx6KLYVIayUWVx6OfPNSz7hcLsaNGwdEo+BtNhtnnXWWafbAABagZCfNZaPpWOunfXl8gls/zx7YxD+r98eWXZqNuYbAqmo8tX8jH9YcBOCq4dO4bEhi89yMcS9ib8tLvW4fmXpJAq3pwHfkn0Q8NbHliKe643NLGYrFEZuWl0gBSlpcTisOm4XzLixGVRUy3Q5s1sT2zzce7ZpAzquHaA3rZNldMfEBeK9qX8IFaJz7qxQ65xDQW1AUULGhKhp1gV1k24vJtI9IqD3t6MHWPrcbJ9g+0JACFA5CKIC+5S3UtCwYPT1pumOapjAo17zSN2lWBwG9a1fGpmooQKY9haZjWQTynObEOKVaC7oNNKfbzJ01tOdOIFD1WccKRQPVCnoARbNjyxoVt2v3N32Jqqo0NDQghODvf/87QJfo+N6IR56iAS1AovoQHI0mJOO9DzEAZfI8lEsSW+2jtKqVNz46jNcfRgiw2zRyLQYWzVwh/MGUi/nVzvdoDHrJtrs4u2AkrXXRCh7fmzSf/1e2E5uqceXwxOX0TnYc+ZOxugsxIqFoy8yWhqJZ0f2NqPYMVGv8BqFPlL4EopHwiqKQlpaGEAKLJSoBQghCoVC3rKXtxCtP0cAWoEPbouLTed3+LZBgAXrjo8N4fB1+GcGQToMnQH62ufl/Muwp/Ghm11QY64gKUEFKOt8uPtcMs5IezZnF8Z1lS2piXAJ6S18CUT+fbdu2AcTSK7ejKArjx49n+PDhPR4br/QlA1qAlEEjgb1dVxYkfuzAMHp+60j6JmIEUVBRFAWBgabY0EUITUmO6iG6v4lQ4wEUqwthhNEDzRiBZjRbGs7hF6AmOPWvy+WK3qteWjlm5AQa2AI0cipkb4OgFwrmomQNQplxccLtuPjsYaz/4DChSLQfbtEUstzm+IucCYQNPx8c/Qlt4You61WsGITJsI3knNz/wKqa14I0wn48+1+LhmEcv83fiL7/NdwJTt3hcDiYNWsWBw4cwO/3k56ejhCCYDDIkCFD+uy6xYsBLUAAOFPBmYq28GrTTBg9NJPbVmYihEAIgaqqrFvX/yTm8SJs6Nz32f+jMejFZbUzK2copS1HAbjpg2dJsViZmDGYa8fMSWil1DLv+93EB8Ag2o1tDh3iiOddxrgXJcymbrYEW3oUn9j2QEsCrekgLS2N6dOnm3LtnpAClCQIIWj1hkixW+hUPARfIIxhiLgVJexr5mRLWogGR7RV5gkHea96P507qL5ImE31pRwtq6DY27MAxWPmpD9dLLO7YZozC9WWihHq2SHSkmZeWE0yIQUoCTAMwd/f2c+RylYcdo2vXTIWAK8/zJMvbEMIOHd6YVxiwfqaOQmo3admm7K6V04IKqLHadx4zZwMdZ1Pjf9zjga2IxAoaFgUG6piAwxyHOMZlnphXK7dXxTNRurYKwi3lqPaUhF6ECPkRfc1ojkzTPWITiYGlAAJw4BD2zAAKvdD1UFocYLFiv7/noDBo1Enz0OxJvbtWXk0mvEQIBDU+WxXLQCtnhDi2FDQx59XMntSAWocasT3NnMy0dPIL0vejy07FY3m4wTIoVpYMf1chvYQjhGvmRNVsXBW7vfjcu7TiWpNLq/n+vp6GhsbycjIoLW1FV3XSUuL+nCFw2HcbjeZmYkNqxlQAmS8+gTsP64SpGU0hPyw7wDs24Sx6yPUlfeiqInzwXHau3ZfUhwWgtBFbJx2S1zEpy+GpWaxesrl7GqpZaw7l1SrnYOtdZR7mxmRmo3DYiXLnoIjwZHdZxJ6sI1IWxVGJIQRasWSkoMta3TCi1+WlZVx4MABACoquo+ftTNixAhGjEjcTPCAESChR7qLT0/UlUFbI6QnLt1ETqaTi88exo59dWSmOzhnWiGvVnxGdroDJT0N3RDMm1WUMHs6k2K1Mzunw7t4bHo+Y9PzTbHlTMMIefCU/G+Xwehwwz7CLWWkjkpsrFp1dfWJdwKOHj0qBSgeKJoFsguhobLvHR0uSEm8P8SUsblMGZvbZZ3FonL1gnEJt0Vyeoh4anucCYu0ViCEQFES16JNS0vrkpCsN453UIw3A0aAANSv3YHY/AYi6IPKA9BaHy1MqKjR6fjsItRL/iXhY0DtlFa1oBuCEYXRZOWBkM7hyhaGD3Yn9Md6PPVBLx8dPUSlr4Xx7nzS7Q4MAVn2FEammpufKGIEOdj2JqqigRDUBXeTYR1JijULt7WILLt5YzBaSnY0Fkzox63PTfjfs7i4GFVVaWlpwel04vf7MQwjVhM+EomQkZHRY+HQeDKgBEhJzUC5YFnXlcfyAWlX32mCRR28+2kZW/dEfWyKR2TR1BrA6w/z8tv7mTQmh0vPGW6KXZW+Zn6x55/ox1LR7/fUd9n+1cKJzC8w5yHXRZh/VN9FwGjqsr4uuDP2eVLGNYxKM6fmm+bIIHXMZYSaSzHCfkTYi+bMwlEwLeG2qKpKcXFyVf2FASZAycyeQx1T1nsPN2L1RzptazBNgLY3V8fEpyc2N1bETYB8Ph9er7fX2bSIEaA5PBToPQJ+o3KQElv34xsaGmJv/3hiceVhcSW+5vqZQnLknZCQl+WMfc7JdGK1dvxpcjPNCykodPZdu2qwCeNl7aiKlRPVW7co3f2WEokR9hNuq8KIBBF6KPo5dOKxmIHCl64FdKKcKBYjzPiGPaSF2zjqyKE1nIYqDD56+n8QKAhFoc6REx0b6oF4ePYCLLpgFJt2VKMbgtmTClj/2l7avCFGTMhnziRzkqsDTMkczMrhM3m1chfeSIgMm5PClAx0YTDI6eaiOHa/UlJScDgcvUZ3AzQGD7C96SkEAl2ECEQasKlppFrzybKPZUzaIjS1+5jea6+9hhpnVwvd34xn/6vHyvI4URQVEfaCaiV1zGVYUhKf2L+5uZmGhgbS09Px+Xz4fD7cbjdWq5XMzMxYeo5E8aUToL48ezUjwqUV/8AadUUkN9SMRY3uN+5oRxeo1FXIluwp3Y6Pl2cvRP185s0aEltWVYX0NDsXzh7Sx1GJYVb2EGZlm29HT2TZR3NhwQNmm9EjoaZDHbNgEX9HR9YIE2o8mHABqq6uZs+ePd3WV1VFk/u7XC5mzZoly/J8UXrz7FWq9mOteKvLunFGd1EZ6qum4Nq7urWC4uXZGw7rfLKtin2lTditGoV5qbR5QxhCsO9II2OHm1eUUHLqaI7eu6eaPfFd18rKvl1QvF4vHo+H9PTElYz+UgpQbwh3DoITjRqAyCrotQt2qvTVNaz15+CPdIzz1DX5sfii9Zv+3z8PkbPxU1Ktvj7PH6+uoeTUsWWNxogE0b21WNxDUFSNcHMpWko2tpzEz0ilpaXR2tp7TmpN02JlmxPFgBIgUjMJL/gO1nefgnAQkeJGZBeiBH0IVzpC1VDsKUSmX3raL91X1zAY6T5Qalg7Zk4CERspWu9lZuLVNSz1NvHEvo8IHCsprKF0mxFTgAnpBXxr1FxUE32VEsmJZud6Zn+nzxXH/nUnnrNzY8aMwTAMmpubcTqdBAIBQqEQTqeT1NRUBg8ejN2e2EH7gSVAgBgxmdCIR0y5dm9dw/c+q2Lv4eYu6wxr1CtaAeadPYmhBb17qH6RrmFfD9MRTyODxImTlfto5cVdVaQdl+84UVPdALoR5tP6X+GJVGNRnOgEKXBMY2LG8oQXGRDH7pnS/llRogPQhnEsxtAcoVZVlfHjx5ty7d6IuwC98MILrFu3DqvVykMPPcSQIdHBzNbWVlatWkU4HEYIwQ9/+EMmTpwYb3OSkgtmDqIgx0nJkRZSHBoFOSm4U6x4AzrZGQ7yO03RJxLRh//P8Ri9pPlMFBvrf0FdcFeXdQc9b2DT0hjrPr2ZB/uanQvW7cFfsRFinf1j90Wzgh5Gc2bhGn0ZqqXnlkYiZueSibgKUHNzMy+++CLPPfccu3fvZs2aNaxduxYAm83GI488Qn5+PgcPHuQnP/kJf/rTn+JpTtKiKArFwzMpHp74CqN9PUwb60t5vnTrCc+RaXNy04SLsGtdf06JfJg8kZoe1zeHjiTk+u34qzYTE53OAq5HszXq/kbCTYex5yafV7IZxFWAtm/fzpw5c7BYLEyZMoXDhzsK3TkcjtiAl81mS9jUn1J75NiPQUGpr0DxtYIewRg1LRoTpkcQg0ad9kHoE/HaB2WU13ii700FCvNcjB7iZvPuOgQKsyfkMG5EYgVqbs4wZmQV4Q0HQVFIs9rxhIOE9Ai+SIh0mxNFUUi3mdNC68xQ1/mUtP79uLUKw1wXJtQO1WLHCIX73EfppfUTT3Rd58iRIwghKCwspK2tjUjk2NiephEMBsnIyEh4Yvq4ClBLS0uXKb2esvELIfjpT3/Kt7/97ZM6986dO3tc7/F4SEnp2XNY++x1LJ/1Ml6y493YR33cWUTmr+z1/Js39yOtx0nY9fG2GspqOgaZhYCKWi8VtR0es+9+Vk12ppOcjJ5nKeJhF4BV1ciwd2w/WbE5Vbv6Y1tnitOvIt06nObQIdIshbRFKhmcMod0W+/+S/G4Zykj5uOv2IgwdBACI9SGotnRnOkYIR/W9CKsGcNPeP7TaZdhGHz00UcxwSkrK+v1HGPGjIkNk5wuu/oirgLkdrspKSmJLffUHF+9ejVz5szhrLPOOqlzT5o0qccR+9LS0l6rPGp7Pu7XudWSjTBvOfTQKktNTWXmzJknZeuJ7Gpu6z15eWe8vnCvAhQPu04Hp2oXnLxtg1KmMyil/wnX43HPLCk5pI39YsnwT9WuPXv20NDQ0G1CIRKJ0NbW1q9zHD16tMdWUPuEQk92BYPBXhsEJyKuHfSpU6eyadMmdF1n165d3UL9n3jiCTRN41//9V/jaUYMkdG/RFrCndOj+MSLWRNye5wXcTk73g8ZaTYG5yVmRqkntjdV8czhzWxqKKMu4OH3Bzbw5P6PqfHLWufJzsmMwyV6ADyuLaCMjAyWLFnCypUrsVgsrF69mnXr1lFUVMSQIUNYu3YtM2fO5LrrriMvL4+f//zn8TSH8MXXY/nsdQh4IBREaaxGCflAgDFoJDjTQFGJzExs+oa8LCfXLhrDwcoW/EEdm1VlzJAMnHaNqrqoA+Lg3BQ0k0o1v1dzgFcqo2+4zxrLu2wr2f0OP5p8KRk2c6u4JhP+6q0Eaz7vsk6xu3EXL4lrKta+JhRqamrYt28fQgjy8vLw+/0Eg8FYXiLDMMjIyGDMmDFYrd1T7MZrQiHu0/ArVqxgxYoVseXOraCe4lLiijONyPnfSOw1+4krxcqUMd1jg4b04f+TKLY39+7CL4C9LUc5K3d4wuxJdo4XHwARbCVYtwdH/qTEGwQUFBRQUGBeUHNvDDhHxGRmz6EmNmyvJRg+Fiyb6SAQ0mnzRmdV0tNsLF8wKuHZ9Eal5XLY29TH9sRnRfyg9sc0hg5iURycl3cPjaH9lLT8HVWxENLb0AmTZhnM2bn/jtOSaPs6+f90WT0wPMVPBilASYI/GOH9LdV0niisawp02aelLcShilZGDUlcsCDA5YPH49Cs7G6uZrQ7l5GpWbxSvhNdCK4aMplcR1pC7TnieZfGULTCQ0T4+bT+cXz6UY5/6Nsilexu+Rszs29KqH32gmkEa7r6Tyn2dOy5yeWFnAwMTAGKhFFLNiLc2Ygh48EwUKoPoNQcRuQUIoYlvpksBPTHmVg3Eu9xrCoKFxWM6ZL7586J5lXG0EWky7JAp8cWByCOy8ecCJyDpuEcNC3h1z0RQghaWloIh8M4HA6cTidtbW2oqooQArfb/eUahE5KIiFsT/8IJRD1sYmMPwfF14JW2uHGr4+aQeSSbybUrBSHhbmT89i442hsndtlJRjWCYaiXTKnXWN0gls/ycjI1Is44nkbT6QGBY0ZWTfSECyhpPUVNMVGRAQAgV1NZ1z6VWabmxQIIdi+fTsNDQ2xdZqmoesdAu12u5kxY0ZCRehLJ0AnjFQOeFEjg2LfXByoRxFGtEDhMURpG6KH4+MdXDm9OIep47JRFYWIYWA59kMI6zqqoqANoBihvlAUlYsG/QzDiKCq0T9kjqOYse4rokGf7cGgCQ5ChWO1wA68gRGMuicomj1amlkY2HOLsZuQhgPA7/d3ER+gi/hAND6zra1N5gOKKxZr15xAqoYwiIpQp3WJJhTWaW4L4Uqx0NwaRFUUVDX6lnI5NEDB6Uj8A7WlsYJnDm/G6KGLk2qxEzYiDHKm861Rc0mzJjbEoF182mkXHDOEpx3vkfdj4gMg9CC6P5rbyV/+CRZXAZozI+F2Wa3Wbi2e41FVVabj+KL0J4+wuvN9tK1vgt1F+PLvoPja0D54AbW5FpGaQfiymyCjeyWDePlClNd6eO2DshOOAZ01JY9p4xKXxlMXBs8e6Vl8ADyR6IN1xNvI2zX7uGrI5ITZlqwIPfiFtscLq9XKtGnTOHz4MH6/n7S0NNLS0mhubiYSiWC1WiksLJQJyRKBMWkexqR5sWWRlk3k6/GtC9ZX17CuyY8WOXHIweYNu6g62L0LGM+uYX+HvHuK8xuIOApn4zv4Nl3unKKCMLBmDEMzsURPeno606ZN67Iu0YUIj2dAClCyoakKfcdPR1HUxPqRaIrKsqHTeL50a49ClKJZCR7rgl00aGxCbQvqbTQESmgKH8GKg/rQXqxKClMzr8dmMS9kxeYuwjL5GvRgC4pqQdFsqBY7Qg+hWqW3+PFIAUoQfXUNvf4w6z8qp6k12jyP6B2Pu6YqOO0arhQr50zNJz+7+484nnl35uQMY06OuW/J42kNV/B+7YPoont3pqp6IwsK1uKwZiTesGOoFhuqJbfLOkWVj1pPDKy7YhgoZbugtQHF7kSk50JtGWrpDozRMyGnEOHKAFdip7pdTitfu3hkQq95ptDY2Nit26oXfI6R21V8glVRvyT74Fre3PFrLDUzTnjenJzTP54WathPuKUMzZWLYnURbjqEEfYhDB3N6sSWOx7bCdJxxIP2XNDBYJBQKBQblDYMg9TUVNLSEutM2s7AESDDwPrSI6gNvcQ1Ve0DQFishC+/GVGY2C5FOGKwfV8DNqvKpNHRxPVb99ZTU+9DUSDT7WD2xFzTAlLNIDc3t8f1IpzG8SNmgaponJN9cC1qIPuELcKcnJxez3+qRDw1+Mo+BCDc0j3nTiTYQsRTg1Z8JZozcaWWDMNg69attLS09LrPkCFDGDMmfkUme2PACJDSVN27+HTeLxJG2/UBkQQKkK4bPPvafvzB6BTp/rIWhOgailFa7WVfaQvXLh6T0OoTe1pq2d1cw4G2ejyRYHS8yjBwqBZUVPwizNSMwSwdNu20X7u3MkNCGGyrf56DzW8RMnzQSY4KnNOYf/G/n3Zb+oMe7F9qEj3YllAB8ng8fYoPQHl5OaNHj054nOGAESDhTEMoCko/ZmtEgoMrm1qDMfEBONoY6HE/XyBCmzdMemr3UsNfhJ66OQBH7QbbMrqGPWTUBbEBjVkd/iIf1x/hSFkpE9u6pnGIVzdHUVSm5V7DtNxrYuvW7VoHwPyhV5/26/VET/fMqhpMzlZxWAwihkLEUKhtif5dh2VHH2xfWGPjJzvQxa5u52w/7+m+Z3a7HUVR+pypdDgcCRcfGEACRIqb8MXfwvLhiygBD2gWhD0FxdtCewUDo2AEImcI+uzefYjiQWqKFUXpiAWzW1UEEAp37WioSjRk43TSVzekydbdaS2zMTr20pzV1WGtwSa6dXvi0c1JBnr7TjoqOxpzSbGECegWdKHQUF8PgF91owvwhm0IRaW3Cc0ves96e5moqorFYomJjMcTTQGcmpqKYRi0tLT0Wd4pXi+TgSNAgBg1jfCoaWab0Q2H3cLieUP5aGstFovCRXOLUBR49f0y2rzRdK0pTgsXn1WE1XJ6x4D6qqa6o7GS/971z375Ap0zZCxXz591+gxLYk6mAu26ddGW2VcWx79ldiLhas8JDcRStLZ7PiuK0mcLKF4vkwElQMlMYV4q31jQNfnYistH97J3YpicVcgdUy5mT1M1Jc21NAd92FQvAsFgpxurptEaCnJW/giWDJ9qqq3JRiTopeXgu4Q9R1E0K8HmChTNgr/hEFZXDim5p3+M8VSE8eqrE9Nl7Q0pQEnCwfJWtu9vIBIxGFnkxjAEaS4bWW4bG3fWYbeqnD9zEE57Yv9kY9PzGJve4b27rvLYD3fW4oTacTz+SBMHm9/GqrmwKE7aQlXoIsK75T8mIkIEIy0Mds1get6/mBIbVr/trxhhP8IQCCNCw65XYh7RAAiDlDxZG2xgCZAQKE01CFVFrdiH8DZB0I+SGq23ZQwaDYPi54/TW//cE3bRGOpo3ja0HEUNR8cODGsO7aGzR6qaKErpXlM8Xv3zZMUQOv8ov5+2UFVsnV+PtsBqfNtj60qaXyWgt3LO4O8m1j49jBH2d9/QKeA51FZjigB5vV5aW1sRQmAYBpFIhHA4jK7rpKYmPv3vwBEgYWB58w9oh7d3jYbvvAvRmmB6LzXBvgh99Z9D4vhaWwpqOJoXqL1GPIAhLKBoqErXUZlEDvYaQhwbMBdEhIHVhMwBQb2ti/j0Rb1/b5yt6Y6qWbE4M4n4O6WxVTQUzYKIBAEFR9aIhNtVUVHB559/DkRfWgBvv/12bGxo+PDhTJqU2GR8X0oB6qmlkRZq4+Ka6Nuxs/iUqFF/jHFGIwqglmzg//w9VyD9Ii2NvvrnB8qa+N93D/a4TVUVjGNZEDPddr5+lTkJtgwh+MmW17C11gFw84fPATA5azA3jz8fSwKFyKG5yXaMpiFw4IT7Dk41Z2A8d9oy2so3oR7ehWKxkjf9q6BqBJvLsaZkY0tLfEbJffv2dVvXeWD6yJEjFBcXY7EkTha+dALUW0sgYnWgo6Id50O7t5MAAeiK1qsXbbxaGqOHZrLssnFsK6kjGIowflQ2mz4+iKapLF48iQ3bqrHbNc6ZVnjar91fNteVUulr4fj39o7GKrbWlzM7b3jCbFEUlflDfkRZ68fYNBcWJYW3P/8Egwij0i/BEGG84TqK0uYyLvPyhNnVxUZVwz3sLCyboy01y7EcQJb8CabYA9EZL5/P1+t2q9UqU7J+UfpqaYhDEzA+ewMaqsDvoVuyCZsD2zfu4uq8ofE1sgcK89MozO+Ix9m9NerU5061c+m5wxNuz/HYtd5/KmZ0w6yqk1EZF8WWXdZoa2hOgbmzOgBGJIi3ejt6yI8eaEUIA1/9wWOD0hE0WwrO7FFxrRHWE7Nnz+bjjz8mEAjEfIKGDh1KKBQiEokwbtw4KUDxRBk5FW3kcdPFx6Yjtav/0wSLOnj7kyPsOdyIrgvsNhW9yYMhBL/4/z4jK93BssvG4bB3LxiXKKZkFzE5azBeWhGATdGwWyzMzBnK1Owi0+xKNoShU7ftb+iBZgD0UPQl11yyvst+gayDZI1PbOvMZrNx4YUXAh3T8FOmTEmoDcczoAQoWdleUsf2ffWxZX9Ax3LMLVoIaGgOsO7tA1yzyNyyLrdNvJB1JdGu6t3nmd/SSEb0YFtMfPoi0HgYIYQp4Q/JxMAJrU5iGlt6mLI9Do8vlABLJF8U1eZCtZw4rak1LX/Aiw8MwBaQiIShrTEamGo79kMJB9GPlkX9hBRQcocm9McxtTiPbSV1fdb8mjnBvDpcyUbY8POPsvtpDVWQZhuM21ZIYyCCoqi8fuQHpNkG4w3VYrekMzv/O7isifORUjUrOVOX0lr6CUbIh6JFZw1dRVNAD2PoISx2N67Bic2fLYTA5/OhqiqRSCSWnD4UCmGz2QgEogHQMid0HBFtjRh//Rm0Rrs7AsAxASIhePrt2DoxaBTqsrtREjQgl+l28J2lUzhQ1oSuC9ypdjZ8eIhQWGf8tMEML0wnP9u8NKOdaQsH8OthdjdVMyFzkCk2fFb7B5qChwBoDh6hOXiEiJgKApqDpTQHS6M7BuGz2t9xQVFix/csDjdZ4xYAYN0fHWtJH3ZWQm3ojBCCzZs3U1NTE1vX3NwMRP2ABg8eTEVF1MF14sSJjBiROB+lAdUFEzs+iIlPjEgPXZvqg1BzKDFGHcPpsDJ5bB7TxuczckgGDrsFd6qduVMGJ434bKorpdrXSnPQz+M736PM02iKHcFI37ltuuyrt8XRkp4Rho63di+tpRsQehiEINhcQai1mmBLJUKPnPgkp5GWlpYu4tMZwzBi4gM9+wrFkwHVAsLl7v++TnNSVG7aWcNnO2uINHlQVYVte4+yv6yJ+iY/E0fncP5M82acSts6BMdAUO5pYmhq4hJrtTMl9xpqS3diEEFBw6pGPcmVYy6mds1NUG9FU+xMyl6aUNuEoXN063PogahIhr0CFDUaC3YMa2oeOZOuQunDteF0YrPZ+swHpKoqhhH1j5N1weKIMvkCaKpFHNkJ4SBkDYY2F3iawWIHRYBmRTnnKpTMxI+5NDT7+WBz9G1kQaAbgnc2lsW8lTbtrGFEUTpF+eaI47TsIg6yEQG4LDaKMwpMsSPLMYKvjf4TzaEyMmxDURSVl7e/jKJYuGzkTTgtmQT0FqyKA6uW2EoUYW9dTHxiCIPO/vdhz1FCnlrs6YlxLE1JSWHatGkcOnSIUCiEEB25m/Lz8xk+fDgHDhxAURQmTEiso+TAEiBVRblwedeV69aBOwft6vjWBesPRg+D0Mev0fUT1w+LF6PTcxmWlkVAj/CdGZeTZTeva2jRHOQ4O1JaqErUR6p9wDnFkviWGYBmSyUqNn1kUVJUNHtiXyKFhYUUFnYIXrsf0OzZs4ET5xKKFwNqDCjZyc1KYfr4roXr5kwuINMdnZkYPzKLoYNOohsZB2yqBbfVYar4JDOaPZXM8QtRbamgWlEtDiyODOxZI7G5B+PIHkXW+IVYHOb+HZOFAdUCOhOYP2co82YV8crfo1UVzptRxHkzitANAy3BbvKSU8OZNRxn1nAALNXRlkZ2gr2ezxTi/ot+4YUXWL58Oddddx3l5eVdtm3fvp3ly5ezbNky3n333XibAoDwtyECPkQ4hNFUC3oE9AgiEkZ4mhBhc2p3t+Pzh2ltC2IY0TEgIQT+QIRw2LyulyEEDQEPdX5PLA2HkSSlmIUwaA6WoYswhohgGAaNgcMEws2m2RNqO0ok5EcPerpt10NejJ5mXhOAruvU1dXR1NQUywMUCoVoaWnB6/USCiXerri2gJqbm3nxxRd57rnn2L17N2vWrGHt2rWx7T/96U9Zu3YtqamprFy5knnz5qFp8QvQMz59FfHhumhmOlUDPQyWaNpT479XRZftTtQlt6MUJr5G0kdbK9m4vRoAiy/64/3juh20eEIoClx89jAmj0lsXz2kR3hsxz843NYAwMjWaCzYI9ve5I7JF2FL0ExOT0SMIK8f+QGecA0tgWiM398OXIcuog/SjNx/ZVzWooTZY0RC1G5+ChHpqGoSblOxuKLjUi2H3sdbvQNFtZBZfDmOzMQFPbe2tvLBBx/EZsLay/S8+eabXfYrLi5m9OjEpQKO669n+/btzJkzB4vFwpQpUzh8+HBsWzAYRNd18vOjs03Dhw/nyJEjjBo1Ki627Nm9iz0bdsQEBwAL1CvRKdxXGBa9Gzrw6nomzNMpLo5/xrq9e/eye/duhIDKo22xP4iiR9MmeI9uja375zu7sepzEmpXWzgAvtZYGg6HP+rD4t96mGdLXsBt7fCcnTBhQtxta7cLIKi30hLKB/KJtEWz+TVu6oiXe4/t7HIGE2bXjq0b0YNdw2qafTqKv451L71IqK322NowyoHXmTr7goT9LTdu3NilhdP+ubq6usu+NTU1hEKhhM2GxVWAWlpaSE/vKHPc2Q+hubm5SzlYt9t9wuJpndm5c+dJ2VJ39GiPdcFS6O4UFtYNjhw5gtfrPalrnAp1dXV4PNGuTZe5E6WHyHchEm5XSHQtzRO2dvTaQ4EgnmDXhFbxtq3dLgCdIBxrMKu2HroPIlp+JlF2BUPhbg+UwwqGouLx+OhczU3XE/u3bA+9aKevnkZpaSl+/4njE08HcRUgt9tNSUlJbLlzrpH09PRYaRCIlgnpLFYnYtKkSSftNCUqRmG8/wKoKqgWOHqsfK7dAalZEApARi7qJdejpCR+lqK0qoV/bCjD6w9jsag4bBpDCtwcbfShqQpfmTuU3KzE+rUAvFq6g3eq92EYBm6rA03VmJ4zhCuGJTaeqSc+q/0Dh1reRQgdpyULpyWTpuBhrKqTCwvvIdOZ2NSnDXteJdhUjqJqqDYX9rQC0kfNQ1GjFTE85ZtQrSlkjL4woVPx4XCYjz/+uMszB1Eh0nUdVVVxuVxMnjyZrKyTc2EIBoMn3SBoJ64CNHXqVH7zm9+g6zp79+5l2LBhsW0OhwNN0zh69CipqamUlpZ22R4PlKKxaNf8MK7X+CIMG5zOt642/6E+nkXDJrMoCcSmJ2bl38Cs/BvMNiNG9vjex5yc2SNxZsev6EFfWK1WLrjgAlOu3RdxFaCMjAyWLFnCypUrsVgsrF69mnXr1lFUVMScOXO4++67+e53v4sQgltvvTWhuWglEon5KKKvgtFJSHtz71S6YBKJ5PTzRZ5J6dkmkUhMQwqQRCIxDSlAEonENKQASSQS0zjjpp3ax8zNiFuRSCTdaX8WT2U+64wToHA4DCQ+daREIumbcDh80kntz7hpeMMw8Hq9WK1WWdZEIkkChBCEw2FcLtdJV1Y94wRIIpF8eZCD0BKJxDSkAEkkEtOQAiSRSExDCpBEIjENKUASicQ0pABJJBLTkAIkkUhMQwqQRCIxDSlAEonENKQASSQS05ACJJFITEMKkKRXamtrufPOO3vdXlFRwfr16+Nqw549e/j4449jy2vXrmXLli1xvaYkcchgVMkps3HjRp5//nl+8Ytf9PsYXddPqvz2unXrOHToED/4wQ9OxURJkiMFSNIrFRUV3HHHHSxfvpwPPviAxsZGKisrufXWW7nqqqtYvnw5Bw8eZPDgwdxwww1cdNFF3H///Rw6dAiAe++9l2nTpnH33XfjcDjYsWMHixYtYujQoTz55JOEQiEGDx7MmjVrcLlc1NbWcu+991JTU4PFYuHxxx9n5cqVhEIhcnNz+c///E/+/ve/s3DhQubNm8c///lP1qxZgxCCBQsWsGrVKgDOPfdcLr/8cj755BOGDh3Kr3/965MSPUkCERJJL5SXl4ulS5eKl156SSxevFj4fD5x9OhRMX/+fCGEEBs2bBDf+973Yvs/+uij4q233hJCCFFdXS2uuuoqIYQQd911l7jjjjuEYRhCCCGam5tjx/z2t78Vf/7zn4UQQqxatUq89NJLQggh/H6/8Pv94qWXXhKPPvpobP+77rpL/POf/xR+v1/Mnz9fVFVViVAoJJYtWyY2bdokhBBi7Nix4tNPPxVCCHHTTTeJDz/8MC73R/LFOeMyIkrM4eyzz8bpdOJ0OjEMI5aZsjMff/wx77//Po8//jgAzc3NRCLR2vELFiyIJZCrrq7m9ttvp6GhAb/fzznnnAPA1q1bWbt2LcAJM+sdPnyYUaNGMWjQIAAWLlzIli1bmDVrFm63m9mzZwMwfvx4KisrT8MdkMQDKUCSfmGz2WKfVVXFMIxu+wgh+N3vfkd+fn63bZ0F5Sc/+QmrVq1i7ty5rF+/nvfeey+utuq6flrPLzl9yFkwySnjcrnwer2x5bPPPptnn302trx3794ej/N4POTm5mIYBq+88kps/YwZM3j55ZeBaLXNQCDQ7RrtjBgxgoMHD1JbW0skEmH9+vXMmDHjdH01SYKQAiQ5ZcaNG0cgEODKK6/kf//3f7n11ls5evQoV1xxBQsXLuTFF1/s8bhbbrmFG2+8kaVLl1JUVBRbf8899/D6669zxRVXsGLFCpqbm5k7dy7bt29nyZIlbNiwIbavw+Hgvvvu48Ybb2TJkiWcffbZzJo1K+7fWXJ6kbNgEonENGQLSCKRmIYUIIlEYhpSgCQSiWlIAZJIJKYhBUgikZiGFCCJRGIaUoAkEolpSAGSSCSm8f8DrnD54dUQnuQAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "# Create a pivot table to check for missing interactions\n",
    "pivot_df = df.pivot_table(index=['ligand', 'receptor'], columns='interaction', values='active_in', fill_value=0)\n",
    "\n",
    "# Convert the pivot table back to a dataframe\n",
    "result_df = pivot_df.stack().reset_index()\n",
    "result_df.columns = ['ligand', 'receptor', 'interaction', 'active_in']\n",
    "\n",
    "# Create a title with key numbers\n",
    "title = 'Unique ligands: {}, receptors: {}, interactions: {}'.format(\n",
    "    len(result_df['ligand'].unique()),\n",
    "    len(result_df['receptor'].unique()),\n",
    "    len((result_df['ligand']+'_'+result_df['receptor']).unique())\n",
    "    )\n",
    "\n",
    "# Plot the boxplot\n",
    "sns.set_theme(\n",
    "            style='whitegrid',\n",
    "            font_scale=0.8\n",
    "            )\n",
    "sns.set_palette(palette=sns.color_palette([region_colors_dict[r] for r in interact_with]))\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(4,4))\n",
    "sns.stripplot(result_df,x='interaction',y='active_in',size=4,hue='interaction',ax=ax,legend=False)\n",
    "sns.boxplot(showmeans=False,\n",
    "        meanline=False,\n",
    "        meanprops={'color': 'k', 'ls': '-', 'lw': 2},\n",
    "        medianprops={'visible': True},\n",
    "        whiskerprops={'visible': True},\n",
    "        zorder=10,\n",
    "        x=\"interaction\",\n",
    "        y=\"active_in\",\n",
    "        data=result_df,\n",
    "        showfliers=False,\n",
    "        showbox=True,\n",
    "        showcaps=True,\n",
    "        color='whitesmoke',\n",
    "        width=0.6,\n",
    "        ax=ax)\n",
    "\n",
    "ax.set_ylim(-0.05,1.05)\n",
    "ax.set_xticklabels('')\n",
    "ax.set_title(title)\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "plt.savefig('./plots/receptor_ligand_interaction_analysis/Club_epithelium_as_source_ligrec_pct.pdf')\n",
    "plt.show()\n",
    "\n",
    "plot_df = result_df.copy()\n",
    "plot_df['direction'] = 'club_as_source'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "result_df['n_interfaces_with_club'] = result_df['interaction'].map(n_interfaces)\n",
    "\n",
    "\n",
    "# Calculate p-values for overrepresentation of individual pathways\n",
    "list_of_pvals = []\n",
    "for i, row in result_df.iterrows():\n",
    "    df = result_df[(result_df['ligand'] == row['ligand']) & (result_df['receptor'] == row['receptor'])]\n",
    "    df = df[~(df['interaction'] ==row['interaction'])]\n",
    "\n",
    "    a = row['active_in'] * row['n_interfaces_with_club']\n",
    "    b = row['n_interfaces_with_club'] - a\n",
    "    c = (df['active_in'] * df['n_interfaces_with_club']).sum()\n",
    "    d = df['n_interfaces_with_club'].sum() - c\n",
    "\n",
    "    arr = np.array([[a,b],[c,d]])\n",
    "    arr\n",
    "\n",
    "    stat, pval = fisher_exact(arr,alternative='greater')\n",
    "\n",
    "    list_of_pvals.append(pval)\n",
    "\n",
    "result_df['overrep_test_adj_pval'] = multipletests(list_of_pvals,method='fdr_bh')[1]\n",
    "\n",
    "result_df = result_df.sort_values('overrep_test_adj_pval',ascending=True).reset_index(drop=True)\n",
    "\n",
    "\n",
    "dict_with_final_ligrec_results = {}\n",
    "dict_with_final_ligrec_results['club_as_source'] = result_df\n",
    "\n",
    "# Save the results\n",
    "result_df.to_excel('./source_data/club_as_source_ligand_receptor_activity.xlsx')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Kruskal-Wallis H-test: test statistic = 19.1446, p-value = 0.0018\n",
      "             Category1           Category2  TestStatistic   p-value\n",
      "0               Immune               Tumor       3.285131  0.001019\n",
      "1           Fibroblast  Luminal epithelium       2.999547  0.002704\n",
      "2           Fibroblast               Tumor       2.936903  0.003315\n",
      "3   Luminal epithelium              Immune      -2.741601  0.006114\n",
      "4   Luminal epithelium    Basal epithelium      -2.052516  0.040120\n",
      "5               Muscle               Tumor       1.943810  0.051918\n",
      "6   Luminal epithelium              Muscle      -1.908803  0.056288\n",
      "7     Basal epithelium               Tumor       1.781672  0.074803\n",
      "8           Fibroblast    Basal epithelium       1.556890  0.119497\n",
      "9               Immune    Basal epithelium       1.536623  0.124386\n",
      "10          Fibroblast              Muscle       0.923079  0.355966\n",
      "11              Immune              Muscle       0.875175  0.381479\n",
      "12              Muscle    Basal epithelium       0.709352  0.478106\n",
      "13          Fibroblast              Immune       0.149240  0.881364\n",
      "14  Luminal epithelium               Tumor       0.003685  0.997060\n"
     ]
    }
   ],
   "source": [
    "# Test whether the number of signficant signaling routes differs \n",
    "# 1) Between all the regions (Kruskal-Wallis)\n",
    "# 2) Between individual regions (Wilcoxon)\n",
    "\n",
    "# Perform Kruskal-Wallis H-test\n",
    "interaction_categories = result_df['interaction'].unique()\n",
    "grouped_data = [result_df[result_df['interaction'] == category]['active_in'] for category in interaction_categories]\n",
    "test_statistic, p_value = kruskal(*grouped_data)\n",
    "print(f\"Kruskal-Wallis H-test: test statistic = {test_statistic:.4f}, p-value = {p_value:.4f}\")\n",
    "\n",
    "\n",
    "# Perform pairwise Wilcoxon rank-sums tests\n",
    "interaction_categories = result_df['interaction'].unique()\n",
    "pairwise_tests = list(combinations(interaction_categories, 2))\n",
    "\n",
    "# Collect p-values and test statistics\n",
    "test_results = []\n",
    "for category1, category2 in pairwise_tests:\n",
    "    group1 = result_df[result_df['interaction'] == category1]['active_in']\n",
    "    group2 = result_df[result_df['interaction'] == category2]['active_in']\n",
    "    test_statistic, p_value = ranksums(group1, group2)\n",
    "    test_results.append({'Category1': category1, 'Category2': category2, 'TestStatistic': test_statistic, 'p-value': p_value})\n",
    "\n",
    "# Create a dataframe from the test results\n",
    "test_results = pd.DataFrame(test_results).sort_values('p-value').reset_index(drop=True)\n",
    "print(test_results)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write the results\n",
    "pd.concat(club_as_source_ligrec).to_csv('./plots/receptor_ligand_interaction_analysis/{}_as_source_ligand_receptor_results.csv'.format(source.replace(' ','_')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Club region as target with all other regions as source"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Club as target in all samples\n",
    "source = 'Club epithelium'\n",
    "\n",
    "results_dict = {}\n",
    "for target in interact_with:\n",
    "    dat_dict = load_from_pickle('./data/region_ligrec_analysis/{}_{}_slides_with_ligrec.pkl'.format(source,target))\n",
    "    \n",
    "    ligrec_results_out,ligrec_meta_out = get_ligrec_results(dat_dict,source,target,reverse=True) ### REVERSE SWITCH DICTATES WHETHER NAG IS THE SOURCE OR THE TARGET\n",
    "\n",
    "    summarized_df = summarize_ligrec_results(ligrec_results_out,ligrec_meta_out)\n",
    "    summarized_df['interaction'] = target\n",
    "    summarized_df['active_in'] = summarized_df['active_in']/len(dat_dict) # Define whether to normalize by the number samples with region interfaces\n",
    "\n",
    "    results_dict[target] = summarized_df\n",
    "\n",
    "club_as_target_ligrec = results_dict.copy()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.concat(club_as_target_ligrec)\n",
    "df = df[df['receptor'].isin(marker_genes)]\n",
    "\n",
    "receptors_to_include = df['receptor'].unique()\n",
    "\n",
    "df = df[df['receptor'].isin(receptors_to_include)]\n",
    "receptors = sorted(df['receptor'].unique())\n",
    "\n",
    "df['interaction'] = df['interaction'].astype('category').cat.set_categories(interact_with)\n",
    "\n",
    "\n",
    "# Create a pivot table to check for missing interactions\n",
    "pivot_df = df.pivot_table(index=['ligand', 'receptor'], columns='interaction', values='active_in', fill_value=0)\n",
    "\n",
    "# Convert the pivot table back to a dataframe\n",
    "result_df = pivot_df.stack().reset_index()\n",
    "result_df.columns = ['ligand', 'receptor', 'interaction', 'active_in']\n",
    "\n",
    "title = 'Unique ligands: {}, receptors: {}, interactions: {}'.format(\n",
    "    len(result_df['ligand'].unique()),\n",
    "    len(result_df['receptor'].unique()),\n",
    "    len((result_df['ligand']+'_'+result_df['receptor']).unique())\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASAAAAEZCAYAAAA39vjlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAAsTAAALEwEAmpwYAABa+UlEQVR4nO29Z2Ac1d23fc1sL5JW3ZLl3nu3sWmhBBcgGBIHG5NAbuAODyUkhEDyhAQIAVKc3Bhy84T0BBJaMCFvaKGEblNs4y73ot5X0vbdmfN+WGklWdXGu7NC5/oizZn239nd357yL4oQQiCRSCQGoBptgEQiGbpIAZJIJIYhBUgikRiGFCCJRGIYUoAkEolhSAGSSCSGkdYCVF5ezpe//OUubd/61rf44IMP+jzvuuuuIxKJJM2u008/HYAdO3awbt26pN0H4Mtf/jLl5eUnfF4oFOK2225j2bJlLF++nH379gHwla98heXLl3PJJZdwySWXnGpzPxUbNmygsbExadfXdZ3rr7+eZcuWceGFF/LYY48BoGla4nlccsklzJ8/nz/96U99Xmv9+vVs2bKlz2OS/XoA/vSnPxGLxQCoqanh9ttvT+r9ALZt29bleU2dOpU9e/YAcM0113DJJZdw4YUX8qtf/ar/i4k0pqysTKxatapL2ze/+U2xadMmgyyKs2TJkpTda9WqVaKsrOyEz1u3bp34wx/+IIQQIhQKidbWViGEEFdeeaU4cODASdkSi8VO6ryBcqK2nag9mqaJd999VwghhN/vF8uWLRPHjh3rdtx5553XY/uJkuzXI4QQ55xzjgiFQid83qmivLxcnHPOOYnt9s9ZNBoVq1atEqWlpX2eb06yWCaV008/neXLl7Nx40ZGjhzJr371K0wmE+eeey4vvfQSNpuNBx98kJdeeonhw4cjhOD6669n0aJFnH766bz33nsAPPzww+Tl5bFmzRq2b9/OT37yE4LBICNGjOCnP/0pDoejx/t/8MEHPPnkk/zP//wP9fX1fOtb36KpqYmzzz6bl156iTfeeINjx45xxx13EAqFsFqt3H///YwbN44NGzbwzjvv0NjYSEVFBTfeeCOXXnopmqZx1113sXnzZiZMmJDoydXU1HDLLbcQDAbRdZ0HH3yQcePG9fpsXnzxRV5++WUAbDYbNpvtpJ7xd7/7Xex2Ozt27ODCCy9k/vz5PT6f119/nYceegghBDNmzOC+++7jyJEj3HPPPTQ3N5Odnc1Pf/pT8vLyOPfcc1m+fDlvvfUWOTk5PPjgg2zevJmdO3dy0003kZ2dzd/+9jeeffZZ/vjHPwLxntvll19OeXk5N954I6NHj6a0tJS//OUvA34uqqomeq9Op5PRo0dTV1fHiBEjEsfs2LGDzMzMLm29PZcVK1Zw1llnce6557Jy5Upee+01MjIyePTRR9m4cWO31/Pmm2/yyCOPEA6HmTVrFnfffTeVlZVdXs8rr7zC17/+derq6ohGo9x8881ccMEFADz99NM89thjKIrC+eefT25uLrW1taxatYqxY8dy2223ceutt/L0008TDAb5/ve/z/79+3G5XDzwwAOMGTOGhx9+mJqaGg4ePEhdXR133303Z5xxBps2beK+++5DVVUcDgdPPvnkgD4fL7/8MkuXLk1su91uAGKxWKJn1idJFMdPTX89oIkTJ4oPP/xQCCHE17/+9cSvW/uvwrZt28SqVatEJBIR1dXVYvbs2YlzO/diHnroIfG3v/1NhMNhsXbtWuH1eoUQQvz2t78Vv/vd77rZ1X7upk2bxDe/+U0hhBB33XWXeOyxx4QQQjz22GOJX4VAICDC4bAQQoht27aJG264QQghxLPPPisuuugiEQgERG1tbeL4F154Qdx0001C13Wxa9cuMWnSJFFWViZ+//vfi/Xr1wshhIhEIiIYDAohhLj22mtFdXV1F/u8Xq9YunSp+MEPfiAuueQScddddyVsuPLKK8XFF18sVq5cKR5//PF+34M77rhD3HrrrULX9V6fT11dnTj//PNFTU2NEEKIpqYmIYQQX/va10R5ebkQQogXX3xR3HPPPYn3589//rMQQog//elPifbOPYaqqipx/vnni+bmZtHa2iqWLl0qysrKRFlZmZg6darYt2+fEEKc0HPpTHV1tTjnnHOE3+/v0v6zn/1M/OY3vxnQc3nrrbcSr+e5554TQgjxox/9SDz11FPdXk9DQ4O4+uqrE72Vu+++W7z00kvdXk/n59fa2ipWrFghdF0Xe/bsEV/4whcSPYz2Yzr3gDp/X37zm9+Ie++9VwghxJtvvimuuuoqIUT8s/61r31NxGIxsWPHDnH55ZcLIeLfn40bNwohhGhpaUk8o2uvvbbP5/ClL31JbNu2rUvblVdeKebMmSN+9rOf9fsc07oHpChKn+2ZmZksWLAAgClTplBRUdHluC1btvD5z38ei8VCYWEh8+fP7/N+hw8fZu/evXz1q18FIBqNsnjx4gHZunXrVr7xjW8AsGLFCv7whz8AEIlE+NGPfsTevXtRVZVwOJw4Z/HixTgcDhwOB7quE41G2bp1K8uXL0dRFKZOncrYsWMBmDFjBt/97ndRVZVly5Yxfvx4AH772992s0XTNA4fPswPf/hD7rnnHu666y6efPJJvvrVr7Ju3ToKCwtpbm7mv//7vxk/fjyLFi3q87UtXboURVF6fT7btm1j8eLFFBQUAODxePD5fGzZsoUbbrgBiM+/dO5VXHjhhQBcdNFF/Nd//Ve3e+7cuZPTTz+dzMxMAM4++2y2b9/OzJkzGTt2LBMmTDjh59JOJBLhW9/6Ft/5zndwOp1d9r3yyiv9zv/0xHnnnQfEP4dlZWXd9n/yySfs3bs3MacZCoUoLi5m+vTpXV4PxOd13njjDQAqKiqoq6vjww8/5MILL0z0MDweT5/2bN26leuvvx6IP7s777wzse/ss8/GZDIxderUxHdmzpw5/PznP+fSSy9l+fLlABQWFvb5HCsqKmhsbGTmzJld2h977DH8fj/f/OY32bdvHxMnTuz1GmktQFlZWTQ3N3dpa+/OA1it1kS7qqpomtbtGp1FrLf/24c5om340C4eJ4LoJaTuz3/+M2PGjGHdunU0NTXxpS99KbHvePt1Xe/VzgULFvDXv/6VN954g5tvvpkf/vCHvYpjdnY2Ho+HJUuWAPEvxyuvvALEP1QQf7ZLly5l586d/QqQ3W5PvMaens/rr7/e7RwhBIWFhTz//PM9XrP9dSmK0usPTW90HhKfyHNp584772TRokWJL1o7O3fuxOPxUFJSckL2QMd72dvnUAjBeeedx7333tulvby8vMvr2bRpEzt27ODvf/87VquViy666JQvqHS2tf0z9/Wvf52zzjqLN954g1WrVvHss88mvme98corr3QZfnXG5XKxZMkS3n777T4FKK1XwdxuNw6Hg23btgHxN6usrIxRo0YN6Py5c+fy2muvEY1Gqamp4eOPP07sczgcVFdXE4lEEnNBY8eOpby8nL179wIQCAQ4evTogO41Z84cXnrpJYDEXwCfz0deXh6KovDcc8+d0HVKS0s5dOgQEP+1yc/P54orruDCCy9M2NgTiqKwYMECdu7cCcCHH37I2LFjicViiVWZSCTCO++8k+gxPP744zz++ON92tbb85k9ezYbN26ktrYWAK/XS0ZGBpmZmbz77rtAvLd08ODBxLVefPHFxN958+YB8Q+t3+8H4j2bjRs34vP58Pv9vP32291+aU/0uQA88sgjRKPRRG+1My+//HI3URrIc+mNzq+n/RlVV1cD0NTUlPi/Mz6fj6ysLKxWK9u3b088s9NOO40XXngBn88HxJ/x8ffozJw5cxLP+O233068z71RVlbGlClTuPHGGxk+fHiPth3P8c/L7/dTV1cHdHy+2nvwvZHWPSCAn/zkJ/zoRz/C7/djsVi47777uvQc+mLmzJmcdtppXHTRRZSUlDB79uzEvptvvpkrr7ySYcOGMWbMGCD+y7Bu3TruuusuAoEAQgi+853vDEjwbrrpJr71rW/x5JNPsnjxYlwuFwCrV6/mlltu4Yknnkh00/ti6dKlvPfeeyxfvpwJEyYwadIkIC4iv/vd77BYLHg8Hn75y18CcZeDH//4x4meTTvf+c53uP322wkEAowZM4abbrqJSCTCtddeSzQaRQjBsmXLOPvss4H48HPOnDl92tbb8znrrLP43ve+x7XXXgvArFmzuPfee1m3bh133303P/vZz9A0jWuuuSYxQVxdXc3FF1+Mx+Nh/fr1AFx66aXcfvvt5OTk8Le//Y3rrruO1atXA3D11VdTUlLSzSXhRJ6Lz+fjoYceYsyYMaxcuRKAO+64I9FT7Gn4NZDn0hvHv54f/vCH3HDDDUSjUSwWC/feey9ZWVldzjnzzDP529/+xoUXXsikSZOYPHkyABMnTuTKK69k9erVqKrK5z//eW6++WZWrVrFFVdcweTJk7ntttsS17nyyiv5/ve/z8UXX5yYhO6LP/7xj3zwwQeoqsrcuXOZPHkyNTU13HnnnT0OwyorK2lsbGTGjBmJtmAwyPXXX08kEkEIwQUXXMC5557b90Pqd5boM0Qyl/BDoVBiGfUf//iHuOOOO5Jyn2Rx/fXXJyaqk43RS8cnQiqfy1Ak7XtAg4Xy8nK+/e1vo2kaHo+Hn/70p0abdEL8v//3/4w2IS2RzyW5KELIhGQSicQY0noSWiKRfLaRAiSRSAxj0M0B6bqeWBE7Uf8RiURy6hFCEI1GcblcqOqJ9WkGnQD5/f5EZLdEIkkfJk6cSEZGxgmdM+gEyGKxAPEXO1B/IIlEkjwikQj79u1LfDdPhEEnQO3DLqvVetIR3hKJ5NRzMlMichJaIpEYhhQgiURiGFKAJBKJYUgBkkgkhpFUAYpGo6xevZr58+cn0oN25j//+Q+XX345q1evZvv27ck0RSKRpCFJXQUzm8089NBDPPXUU932aZrG+vXrefzxxxPZ05544olkmiORSNKMpAqQoiiJNJ3Hc+TIEUaPHo3b7cbtdhOLxQiHw3JpXSIZQhjmB9Tc3JzI9wvx/M5er7dbYq3eaM/2J5FIBi+GCVBWVhatra2J7dbW1n4TbXdm+vTpsrckkaQB4XD4pDsEhgnQqFGjOHLkCIFAAL/fj8lkkoIikQwxki5At9xyCzt37sTpdLJ9+3Y8Hg/nn38+Y8eO5aabbuJrX/saiqLwve99L9mmSCSSNGPQZURs7+7JIZgk3RBC0HLkPcJNx7BmFpE19iwU1WS0WUnn03wnB10wqkSSjpSWlrJz22Ziwaa2lkbMH5cTisUDNI8vgAgwderURNWLoYoUIInkVCH0rptCJxCIV8LtSYAkUoAkklPC5MmTmTh+LA27nifqq8XsyCZ3xqU8///Fi0xedtllBluYnkgBkkhOEarZSt7ML6FHg6gWO4oiQy37QwqQRHIKURQFk1UOtwaKFCCJJAX4Kj4h3FyBzTMCd3H3GvdDFSlAEkmS0aNBWo68B0C46QhmWwb23DEGW5UeyEGqRJJkhK512Y6Fmg2yJP2QAiSRJBnV4kC1OOL/W1048sYbbFH6IIdgEkmSUVQTBXOuIBpswuLMQTVLD/52pABJJClAtdixWYqMNiPtkEMwiURiGFKAJBKJYUgBkkgkhiEFSCKRGIYUIIlEYhhSgCQSiWHIZXiJJNkIHV/VDrRgMxZ3Po78CYZFytfV1bFz505UVWXmzJlkZ2cbYkc7UoAkkqQiiPrraTlUk2gJe4+RPfHzKbOgtLSU3bt3A9DU1ISu62iaxt69e8nLy+t2fCozNUoBkgxKhBBsqfsTFa0fkeuYyGnDbsCkWo02qxtCjyH0GKAk2kKNR42zpy0FvKZp6Lrez9HJRwqQZFBS7vuQfU0vAuBvrSPHPoYpOZcYbFV3FMUMigp01H6wZg5LqQ2TJ09O9GiOHj3Kzp07qaqqwu12G56pUQqQZFAS00NdtqN60CBL+kFRsLjycRUXEws1Y3UX4jIwH9CoUaMoKSnh+eefN8yGzkgBkgxKRmYs5lDLf6gN7CLLOoKJnuVGm9Qrimoia8zpRpuRwGRKn1JBUoAkgxKTauW8EXcT1QJYTDIF6mBF+gFJBjVSfAY3UoAkEolhSAGSSCSGIQVIMmgJxppoCVcYbYbkUyAnoSWDkmOtG9lYuR4djTGZn+O0ohuNNklyEsgekGRQsqfhH+jEq00cbnmTQLTBYIs60LUIsXBrwutY0juyByQZlDgsORA+BIBZsafNaliktZqGXf9CaGHsOWONNiftkQIkGZQsKPxvTIqFsNbK1JzLsKgOo00CwFe+FaGFAQg1HkJouSim9ItRSxekAEkGJQ5zNqcX32q0Gd1QLfZOW0pbHJikN6QASSSnkMzRi9FjEbRQM66iGSg1e4w2qU/q6uooLS3FYrEwY8YMXC5XSu+fdAF6+umn2bBhAxaLhfvvv58RI0Yk9j333HM8/vjjqKrKRRddxFVXXZVscySSpKKa7eRMXtqpJb0F6OOPP0bT4pP527ZtY8mSJSm9f1L7h16vl2eeeYbHH3+c73znO6xbt67L/kcffZTHHnuMp556iqeeeopIJJJMcyQSSSeEEAnxAQz5/iW1B7R9+3YWLlyI2Wxm5syZHD58uMv+sWPHEggEALDb7WkVpSsZPBxteY/64AGaI0dxWQqYX3gNJsViiC3RQBONe15EC7fiHj470a7HIjSWvkikpQp79kiyJy1DUY39vCuKwvjx4zlw4ACqqjJp0qSU25BUAWpubiYrKyuxfbxfxNKlS1m5ciUmk4lrrrnmhARo586dp8xOyeClwfI+Dba3u7RVNO5gVPC/DLEnO7IXu+4FwFe+GX+rE6GY2Lf5FTJjca/tUOMR9nz0b4LmAkNsBPD5fAD4/X4KCgpQFIXKykoqKytTakdSBSgzM5O9e/cmtlW1Y8Tn8/l49NFHefnll7FarfzXf/0X559/PsXFxQO69vTp07HZbKfcZsng4o2yf0Gga1vU1Mi8efMMsaextI5QQ1Ni2+lyoahmSkpKaDlyLNE+eswYnAWpybvcE0ePxtPCnornFA6HT7pDkNQ5oFmzZvHRRx+haRq7du1i1KhRHTdWVSwWC06nE6vVit1uT6iyRDJQCp3dswtm2UoMsCRO5qjFmJ1x35+MkYtQ1PhvvHPYdGzZo1BUC/bc8TjyJhhmYzqR1B6Qx+Nh5cqVrF27FrPZzH333ceGDRsoKSlh4cKFXHrppVx++eUoisLs2bOZOHFiMs2RfAaZlnspGZZC6kMHaQ63zwEZM/wCMDuyKJizuqPho3ivRzVZyJ16kUFWpS9JX4Zfs2YNa9asSWx37gVdffXVXH311ck2QfIZZ2TmEkZmpnb5WHJqkG6aEonEMKQASSQSw5ACJJFIDEMKkEQiMQwpQBKJxDCkAEkkEsOQ6TgkkqQi0EIt1G55AhQFULC688kYdRoma3pkcTQSKUASSRKJhVrQI35iwY54kVigHi3iJ3faxYbZ5fP5iEQi/Otf/8Jut+NyuRIBqR6PJ2V2SAGSSJKJHuuxORZqSbEhHZSVlREOhxPboVCIUCgEQEtLC+effz6KoqTEFjkHJJEkEdXac4ZBd/GMFFvSgd/v73VfOBxG1/WU2SJ7QBJJElHNdizuQnKmLATVBAhMVjcWZ7ZhNo0ePZp33323x7JBY8eOTWleLilAks8EEc2PRXWgpGESeEU1Yc8Z1f+BKcJut5OdnU00GmX69Onk5ORgsVjQNA23251SW6QASQY1utB4r/KXlPs+xGnO49wRPyTDWmS0WR0IgRbx0XJ0E67iWZgs6VE+SFEUrFYro0ePNtQOKUCSQUVpaSm7d+9ObEd0P95wGD28gCZgg/NZMqxdk9pNnTqVyZONSf4VCzaix8L4yjcTbjpG/uwvG2JHupJ+/VWJ5ARQ2j7CesSKHrGiKOmVV1zXoon/o/46hEjdBO9gQPaAAKFriPefR9QeRZm0EHXa6UabJOmFyZMnd+vN7Gx4lo2vHMSs2lj9pa9iNaW2tlVfqGY7erSt8ELOGMPnqAKBAKWlpbS2tuJ0OqmoqKC8vJyMjAwmT57cJW1yKpACBIgtryI+fCH+/5FdiJwilCJZ13uwMD33i+yzbwBIK/EBMDs86BY72RMXYM81/jO1ZcsWvF4vkUiEWCzG1q1bgXiBQrPZnPKspHIIBtDc0GlDIFoaej1UIjlRVLMdR/4Ew8vwAIkyWEA3f5/O+1KFFCBAmXEG2NricnKHo4yebqxBkkGNFvbjr95JuDlehkfoMfxVO4n66gy2LO7n047dbk+UzTKZTF3SJacKOQQDlIJRqF+7H5rrIK8ExWI12iTJIEWPhqjb/nf0SLzCixbOQAu30nzoLVBUcqdfgi1zYKWnksH48eMZNmwYL730EiaTiSVLltDS0oLT6TSkzJUUoDYUZwY4M4w2Q3IChGLNHGr+D5qIoQsNta0aakTzUenfittSSJ4jtXMaUX9dQnwA9GgQEIACQifcVGaoAAG43e6Et7PJZCI72zivbClAkkFJWGvlxSPfJqw1A9ASmkOOfRxRPci/j/5fWqNVgMJpw25kTNbZKbPL7MxBMdkQWjzYUzVb0SIdS/HWjMKU2TIYkAIkGZQ0BA8kxAdAoBPV/XhDR9vEJ95a1roxpQJksrrIm3Epwfr9mB0eTDW7US0O3CWjsGYUYc8ZnTJbBgNSgCSDkixbCQpmBB3pLsyKnQxrERbVSVSPr+jk2Mel3DaLKxeLK7dtazeKyUrmqMUpt6Mzra2tfPDBB4TDYRRFoampCUVReP/993G5XEyePFnOAUkkA8Vlyef8kT9iR8NT6HoMxToKs2rHbs7ivBF3c6jlTTIsw5jgWWq0qUDc2TXSUoUW8WN2eFI+FHv//feJRuNDQSFEYgm+sbGRxsZGQqEQixYtSqlNIAVIMojJc0zgnJI7Adjw8YZEe7Z9DPPsY4wyq0cadv+LSHN5YttdMp/MUan7wsdiPSdGa6evHEHJRPoBSSRJRuixLuIDEKjemVIbCgoK+txvVFS87AFJJElGUUwoJitCiyTazI7ULn0vWLCAqqoq6urqcLvdtLa2oigKp512GlarlczMzJTa044UIIkk2SgKudO+QGvZR2ihFizuQjJHn5ZyM4qKiigqiudK+uSTTwDIy8tLuR2dkQIkGZQcaX6X7Q1PYFGdnFH8bQCieoD3Kx/EZSlgeu6XMKnp49FuzSgkd+pFRpuRdkgBkgw6/NE6NlavT2y/duwH6OLzeMPH0Fvj0d2aiDC34GpD7NO1KHoshMnqBiEAmQOoN6QASQYd/mh9l+2w1oIqYohOX/TmcPnxp6WEqL+Bhl3Po0eDWFz5RFprAfAe+A+e8ecYYlNfaJrG/v37CQaDjBo1ipycnJTeX66CSQYdOfax2Ewdk6bFrrmYVRsWNZ5vWUFlnOc8Q2zzVX7SFv8VjwtrJ1Cz29BaYL1RWlrKgQMHqKioSDgqphLZA5L0SygW5bel73K4tYE5eSNYO34haooK1/WEWbXxhbGPcND7Gg5zNiMzl7Dhww1k20azpOQCnOYcMm3DDbFNNdt73qGYUEyW1BozAHy+jsBZTdMIBoMp9YhOugA9/fTTbNiwAYvFwv3338+IESMS+xoaGrjnnnvwer3k5+fzi1/8ItnmSE6C1ypK2dkUj696t/og07OLmZM3op+zkotZtTEp58LjWhWGuYwr+AeQMWIBejRALNCEPXcs6uFPEEIje9Ln06YiRmdGjBhBfX09Qgg8Hk/Kl+OTKkBer5dnnnmGJ554gt27d7Nu3TrWr++YPPzJT37CbbfdxsiRI5NphuRToh2XSD0mE6v3imq2kj3x84lts/MYAI40SMfaE8XFxWRkZBAKhcjJyUl5Tuik3m379u0sXLgQs9nMzJkzOXz4cGKfpmkcOnSIhx9+mCuvvJIXX3wxmaZIPgXnDZ/EaHcOCjA7t4S5Bvd+JKeWjIwM8vPzU1oRtZ2k9oCam5sTKR+BLqVgGxoa2Lt3Lz//+c8ZNmwYV1xxBUuWLMHj8Qzo2jt3ptaVfahzAfkIRx5KUOGTLVuNNqcb7XMZmzdvNtSOjOhRXFo1McVOk2USvtZWLMJHxXv/C8TnzaKKm0brJIRi3BRsujyvpD6BzMxM9u7dm9ju3L3LysqiuLg4kaN22rRpHDt2bMACNH36dEPSB0jSj/LWD9lcWgmAe0wOe5teJKoHmFNwNWOzPnfK7/fOO+9QV9c9v7PTHKUoN+4iYBFB1OZSiIVRTAIFhXhmRLCKVmINpZT7e87AmZ+fz5lnnnnK7e7M0aNHAZg3b96nvlY4HD7pDsEJCdCxY8eoqKhA07RE2xlnnNHr8bNmzeKRRx5B0zRKS0u7JL222WwUFhZSX19PdnY2+/fvp7jY2FSVksHJpur/RWcCAFvrHkcQ/3x+VP1rRmYsxqye2h+quro66uvru/nMHF9lQhdQmGUi09p9zqxzSozONDY2nlJb050BC9A999zD1q1bmTRpUpeeTF8C5PF4WLlyJWvXrsVsNnPfffexYcMGSkpKWLhwIbfffju33HIL0WiUiy++2PC4FMngRCB6aU8eOTk5rFixolt7sHIz4bo9mGwZjJt6LhPMNlr2PIeIBkBRQVExO/OZNuMcZpi7C+NQmwsdsABt2rSJF1544YRnydesWcOaNWsS2517QTNmzOCvf/3rCV1PIjmeRcP+Dy/xJgowK28Ne70vENWCzC24+pT3fvrDUTwPR3HXYU3W9MtTasNgYsACNG3aNCoqKrr48Ugk6cDIjMXkO+J+SlNyL2FK7iUGWyQZKAMWoEOHDnHRRRcxbtw4rFYrQggUReHJJ59Mpn0SyaAl2lJJsOID9GgQRTUhdA30GCZHNs7Rn8Nkk2WgBixADz/8cDLtkEg+Uwih4z/8BuhteZg71m3QAvUEy97HPT498lUbSb8CFIvFMJvN5Ofnp8IeieSzgdAT4tPj7lhqgz7TlX4F6MYbb+TRRx9l2bJlKJ0CENuHYK+//npSDZRIBoImougiiqZH0iIRmaKasRfNJVS1paed2Ivnpt6oNKRfAXr00UcBeOONN/o87pNPPmH27NmnxKhUIyoPIPZ8ALlFKLPO6SK0kvSn0reVxtABBILXy+7ivBH3pIUI2YfNwpo3qc0fQAfFhNDCqGZHyiPjw+EwBw8eRFEUMjMzaW5uRlEUtm3bhqZp+Hw+bDYbs2bNwm7vJaI/CZwyT+h77rmH55577lRd7pRTWlrK7t27u7UHfK3QXI9TRIBjsP0wuLOZOnUqkydPTr2hkhPmYPOrCV+ghtABGkMHyXdOMdiqON3Sc/Tg+3Oq6M1DG8BsNidcaIQQiTI9ZWVlXY579dVXE/XDjicZHtqnTIA6x3kNJgKBAGDC2d4QCRlojeRkcFuKgHhogapYcFqMc2iNNB4k4j2M2ZGHbdgsFEUhWPExkcYDqFY3rvEXoJqS0zvrzUMb6NKrVxSFjIyeV+AURUmph/YpE6B0H7ZMnjy5xx7Nhr//HWqOcEn4ACgK6uduRBk/xwALJSfLzLzV7DD/Dk1EOWv4HbgsxiyYxPx1BI6+Hf+/uQzFZEW1ZRCu3QGAFgviP/gaGRO7e1CfKnrz0N61axc1NTWJbbfb3eP5Ho+HuXO7z08ly0N7yPeAUFUoHIUy4yyU7EKUglH9nyNJK0yqBbclXuq4yDXLMDv0cGuXbS3SgjhuJUyP+DCCqVOnUlBQkOj9HDt2DJvNht1uRwhBa2srNpuN4cNTm0nylAlQT6o5aFBNqJMWGm2FZJBjzhyOastED7eAasGaMx7V4iJUvS3hCGTLn2qIbYqidHGlmTBhQpf9hYWprVXfzoAF6OjRozzwwAPU19fz97//nX379vHWW29x3XXXAfDDH/4waUZKJIMB1WwjY9IX0IKNqLYMVEt8ZjFz2peJeg9jcuZhdkl/us4MOLL0+9//PjfeeGNi9nzChAk8//zzSTNMIhmMKCYLZndhQnwAVIsdW/4UKT49MOAeUCgUYsaMjoTfiqIYksJxKKDrgu376mjxhZk+IY+crPRLZi7pHSEEkcb96KFmLNljMTtz4+1ahHDdboTQseVPQ03iknx/NDY20tDQgMfjMTTKYcACVFBQwJ49exKrXU8//bSMjE8SGz+p5IMd8ejuXQca+Nql07HbZAWlwUK4bjehig/j/9fvJXPKpahWF/7DbxJrrQAg1lJJxiRjSjU3NzcnasOXlZUxc+ZMw3JxDXgIdu+99/K73/2Ouro6zjzzTDZu3Mg999yTTNuGLFX1HSslwXCMZp+MGxpMaJ0KEqJH0ULeeHugo10L1Bu2ctzS0tLndioZ8M9qRkaGrNuVIiaMzOZYVXxJNyfLTk5W6lzjJZ8ei2ckUW+8AoxicWFyxnsXlqxRRBr3t/0/0jDfuZycHEwmE5qmoaqqoZlIByxAF110EdOmTWPp0qV87nOfS2m8yFBj1uQCsrPstPojjBvhwWKWc22DCWv2WFSLCy3cgiWzJDHX4xh5OuasEhA6Fs9ow+xzuVwsWLAAr9dLZmZmr06JqWDAAvTKK6/w8ccf88orr/DLX/6SyZMns2zZsh69LiWfnpFFqa1QKTm1mN2FmN1dfWsURcFqoPB0xul04nQ6+z8wyQx4DkhRFBYsWMCdd97JY489hsVi4dvf/nYybZNIBiVC14g0HiDcsJ9Q/V6i3nicWrS5jHDDPoQWMdjC9GHAPaDq6mr+/e9/8+9//5tAIMDSpUu5+eabk2mbJA2I6hoRLYbLImuwDRT/kf8Qa+4aZW5yDUPzVwMQrttDxqSLUZTUlkFORwYsQDfccAPLli3j/vvv/8zUcheBFmhpAFVFaDEUk1zq7sxebw2P7H6LkBbjvOJJfHncpy9i91lHCEGsubxbe+cVMD3YiB7xYbKldpgdCoWorKzEZrNRXFxMU1MTXq8Xl8tFIBDAarVSXFyc0snxAX/jNmzYkEw7Uo4QOvozP4eW+DhYvPYXlKX/ZbBV6cU/j24npMU931+v3Mt5wyeTa3cZbFWcSv9WttX9FbNqZ2Hh9YQ0L4FoA2+W38eiYTfgMGcbYpeiKJhc+Wj+2i7tqi0TPdQUP8bi6uIpnQo0TWPz5s2Ew3GXjqamJmpra7sdFwgEusWJJZN+Bei73/0uP/nJT7j88st7TMk6aKtihPzQUAnm8QCIigMGG5R+OMwdWftMiootTTzfNRHlvYpfEBPxL9P7letpiWQBUOXfxtbav7Ck+BbD7HONPZ9I/R6EFkMgMFlcWPImEm3Yhx4NYcubiKKmtrcdCoUS4gNxZ8Se6K09WfT7FG699VYAfvnLXybdmJRid0PRWGjrGStjZxprTxqyZtwCovomWiIhVoycjtuSHq4XQmjERMdEbkT3A1mJ7ajuN8CqDlSzDfuw2d3ajYqEB7Db7TidzrYEfJCXl0dlZWU3Z8jc3NyU2tWvABUUFADwi1/8opsI3XrrrYNWmBRFQf3it+GpJ0BVUc5On+qV+4408sYHxzCbVJadOYaSQmPqR+XaXXxrxnmG3LsvzKqdmXmXs73+SUyKhXkFV/O66UOCmher6mZa7ipD7dMjfrRwM2ZnHkpb9kOha8T8tagWJyZ7Vj9XOPWYTCbmzZtHbW0tNpuNvLw8iouLaW5uxu124/f7sVqtKY8L61eAdF0nFotx6NAhotFoQjF9Ph979uxJuoHJRLHawRX/MKRLRkddF7z87hFiWjwt5qvvH+Vrl0432Kr0Y1ruF5ngWYaqmDGrNjKs5bgpZOW4OzCpqU343pmYvw7fgZdBj6FaM3BPuhjFZMF38BU0Xw2g4Bx9NtbsMSm3zWKxdEk4lpGRkUjN6vF4Um4PDECAHnvsMf785z9TW1vLsmXLEgLkdru54oorkm7gUKRzt3jQZppMAVZT1wlxBZOh4gMQaTwAenziXo+0EmutQLVltYkPgCBSv9cQAUpH+hWgq666iquuuoq//vWvrF27NhU2DWlUVeHzS0bx+qb4EOy809LT5SGma7xUtpuGsJ+zho1nbKZx8UTpRNfhlYJqy0S1uEA1J4RJNWAIlq4MeCq+qamJlpYWMjPjvgvNzc089thj3HTTTUkzLhkIXUfseheCPoTJDK2NiWFYujB1XB5Txxn7he61jFHbJGZA1WkMx///u7KNMRm5zJg2fciXMrLmTUFoMbRgAxbPGMxtgajucRcQrtuDanVhLzKm6IGu6xw4cIBAIEBmZiaBQIDm5macTiezZs1KlO1JJQMWoNdee62L2GRlZXVrGwyIN59AfNKpyKJ5PARaELqGoqbHMjNAc2sYm82EgkIwHCPLbU2Leap2AQo7Oj46uhBEda23U4YUiqJgH9Z9RbWn2LBUs23bNpqa4r5IncvshMNhNm3axJIlS1Ju04AFSNM0AoFAIoDN5/OhaYPvQyfKSrs3RsMQaAG3Mc5rnRFC8NI7hyk93IjJpKAoCrGYzrgRHr5wzriUiVCvZYzaHFJHnDGX3+/diEAw0p3NqlkXYEkjAU8nhBZFj4VRFFDM9pT7ALXTl49PKGRMPbwBP4mvfOUrrF27lhUrViCE4MUXX+Sqq65Kpm1JQRk1DdFQ2bXRYgenMcOw44c6sZhOdYMfS1s5F6FaMANH98JTdVuZPWtGWgxzFhSMptjloTHsZ1JWoRSfXoj5qvEdeBVEfP5HMdtwjVuaSNOaSjweT68FBo2KjB+wAH35y19m5syZvP7662RlZbFu3boBGf3000+zYcMGLBYL999/f7c0rq2trZx//vncc889LFu27MRfwQminP1lKBgFIR/CbIVPDoAzC8WA8W9PqKqCggKivZ6Upcu+dGK4y8Nwl8doM1JOIBDA7/cPqFifFmxEaJ1rg4VQ9r+Cye7p8fiGhgZcruSEu8ycOZPDhw/j9/vxeDz4/X68Xi9Op5OZM41xxD2hOaAHH3yQiooKSkpKeOCBB5g0aVKfMWJer5dnnnmGJ554gt27d7Nu3TrWr1/f5Zjf//73zJqVumJyiqKiTF3c0XCgIWX37omehjr7jzbx+qsvoKoKw8YswheIMntyAVPHpf5XU/Ip6Sni3aAoeFVVGTdunCH37o0BC9BDDz3Ek08+ydq1a3n++efZvn07f/7zn/s8Z/v27SxcuBCz2ZxQ387U19dTVlbWpdqGBCaMymZHTrx3+YVzxhtsjeR4nE4ndrt9QMn49GiQYNn7xAINoCiYnfk4Ry5JeEgfz4svvmjIapRRDFiALBZLInVjLBZj5syZ7Nu3r89zmpubycrqmFs53qnu17/+Nddddx2vvvrqidgMwM6dO0/4nJ7w+eIJ4Ddv3nxKrneq8Da3ogvY9OHHWEzpM/RK1+eVSrt8Pt+A50xUiwPX2BMLZ/H5fCf1Ok7ErpPhZO3qiwELUH5+Pi0tLZxzzjlcf/31ZGVlJeLEeiMzM5O9e/cmtjsre1lZGS0tLUyePPmkBGj69OnYbCeWJEsIHfHW0whvLeSNQMnO56jbDYEW5qhelCmnoaRBuol9RxrxheICu6fCxtqLpmK1pMck79Gj8ex+8+alV26gVNp19OhRdF1P2vXdbvdJvY7+7KqpqSESiVBYWEg0GqWhoYGMjAyyswe2+tubXeFw+KQ7BAMWoF//+tcAfPOb3+SDDz7A5/Nx5pln9nnOrFmzeOSRR9A0jdLSUkaNGpXYt2fPHo4dO8Y111zDsWPHcLlcjBs3Lqm5SPRn1kF5myAe2oYAcEyDaBjxn3cQO95GvfKHhvsDlR7uWKloaglT0+BnxDCZI3qwEG0+RuDouwDYi+cTqS9FC3mx5ozDMWKJIf5chw4d4siRI0D8xz8ajSbcaGbMmGFYccKTckhYtGjRgI7zeDysXLmStWvXYjabue+++9iwYQMlJSVccMEFXHDBBQA8/PDDTJgwIfmJkGqOdG+LdcrPW18O/mbIyEmuHf1QkOPkSNv/FrOKJyM90mCkE97wMar9O8i0FuENH6M1UonV5KY2sIcjLW9jM7kZ71mKy5J6j/LAsfcQWjz3TrB8I4h4ryTSsA9L9hgsGcUpt6nz8vvxPj+NjY2DS4BOhDVr1rBmzZrEdudeUDspyy2dVwJVB7u2mSwdIpRTlBZhGYtmFrFvh51YTGf5BRPJcPU8YTlUaQlX8O+j/xdNdCTYCmqzCGpeXi97IdF20PsGF459EJvJmHQmcdJj/s7j8SQKENpsNqLRaGK4ZlQkPKRAgNIJ9fLvor/2F2iug7wSlJxhcKARAi0ok2ahTD/T8OEXxN353c64/09RvnE1m4zmnXfeoa6urlu7ln0Arbj/arFhvYX/79XHUQM9h0Dk5+f3O41wMjhHnkHg2HuAwF68oG0I1oQ1Z7whvR+AcePG4XK5CIfDFBUVEY1Gqa+vx+12D47ChJ8FFFXFdMHVXRsPbgBXFurC9KxvtvtgAw3eIBNHZ1OYa/wEeSqpq6ujvr6enJzjhsS+HNBVUDsmXO3F1d0vELMgghk9Tsz25hF8KrBkjSBrxurEti3XeFcKRVEoKipKbNtsNkMLErYzpARosOEPRHn53bjv1CeltVx1yTQy3akrjxOKRdlw5BMaw37OK57MjqYKjrbEeyS/2vEmmTY7LdEQSwrHMTdvRN8XO0lycnJ69LdpCi+iNrQDl3kYLdEyfGOrKHDMxK56KAu8h1V1MS5jGe4xw3q87kC8mCXJRwpQGhOOaonSkdGYTn1TMKUC9MzhLbxbHZ8z29VYhY5gTJsv1w5vRzzdrqYq7pp7IcOcqVupy7aNI9vW7tV7Wpd9w5yzU2bHQNFCzQTL3kdoEezF87BklhhtUlowJAVIRMPou96Pr3oFdDBb0F58FMx2lMJRKOPnoKTBZLTDZqalLYzI5bAwLD+1Q7CGUEdyd53eMzPqQtAY9qdUgAYbgWPvofnjWRH9h/9D1ow1KY2KF0JQW1uLpmkUFhaiKAo1NXF7CgsLDfO+HpICpD/7S6hsK8PTVpanfVvsBPHxS6hfuRvF6jDIwjgOu5nV50+mwRtk9PAsnPbUphv9XPFE9jXXogmdCRn5HPF3nTcxKQqaEIx25zA+05hl3MFC+7I8EM+MKJLnyNgT+/fvp7w8XjCxqqoKm82WqAtWX19vWDjUkBMgEfJ3iE9vNNdDQ1W8bI/BFBe4KS4wZrJwdm4J986/mJZokJHuHMKxKE+V/R2TonLl3GW4zTYawwFGuLNlOo5+sBfNJXDkLRAatmGzeo0F+zT0FaXf3NycmIyvqqpCUZREaFR1dTVlZWXdzulMsqL0h5wAYXNCbnG8KGFvODLA03eYyamktjHAe1srMKkK/kCUxpYQlmAQi0nl7//ey8RROcycZEwPI9fuSlRDdVpsZLX1Cktccfd9j82YPDLHUxn4iL3N/8QXq0QgyLdNZXHBbUablcDqGYVlxhqE0FDNqXcsNZvNRCJxfzeTydSW6C6W2GcUQ06AFEVBXXU7+scvx+eAWjPAbAEToJqgeDzqtCUojtT0OoQQPPfafvzBaJd2LRwjBDRVtXKsqpUcj92w+mDpji9azccNjyDoGNbUhnewx/ssUzxfNNCyrigmCwrJG0b3FaWv6zrl5eVomsbw4cNRVTUxJCspKelXhJIVpT/kBAhAcWZgOquteF1bPiPTZd8zxBZNF93Epyda/ZF+j0kVQS3KwZZ6NF1juMuDy5K6lbmeCGlNXcSnHV+sB9+gFCKETqh6G1HvURRFxZIzDnvBNENsUVWVkSO7VlgZPXq0IbZ0ZkgKUDphNqnMnVrIlt013fapbUGLuR4HY0uMX5UTQlAZaMYXDfP2tn8D4LE6uGP2BeTYjHOSzLFNINc6iYbI3k6tKhMzv2CYTQDh2l2Eqz9JbGsVDagWB9Zs4+cW0wUpQGnA5xaMYMaEPEwmFVWBilofu7YcBQXO+txUcjx2zCbjk1Q1hP34ol1DILyRIFvryzhvuHF5qlXFzJKC22mJlOOLVhMRPoodC7GbjXUL0MMtA2obygw5AdLLSmH/ZoRqjlfC8GpgcyL2fYwQAmXCXEPiwXI9HUv+mW4bez6J934KctNjkhfAbbYllt47U+Awfm5KVcx4bKPx2EYbbUoCa+4EIo0HQcTTXigmGxaDej+hUIjGxkYyMjIQQuD3+8nJycFms9Ha2kprays5OTnY7amdIB9SAqTveg/xyh+6NprHg68J/V9vxbcnLkC56PrUG9eGPxjlX28epLY2nuHv98/uYPmZYwxbiu+M3WxhuMtDUzjAuMw8nCYrM3OHMyNneP8nD0HMrgIyp30JLdAAiorJmWvIClgoFOKjjz4iGu0612i1WpkwYQK7d+9GCIHFYmHBggUpFaHPnAD1FkENsLDmY4p63NOBtu9j/tFHov1kRVC388H2KipqfZjbehnNvjCvbjzCVZdMT9o9TwS7yUKRM4vLZl1gtCmDAtXiRM0ythfr9Xq7iQ9AJBKhuro64Q8UjUbxer0MG9Zz/Fwy+MwJUK8R1ECjNYuiYNfJ3sl6V+/eJmtWr2ktkxlB3Y6udw956KlNAhHdj1V1oYkoKmZ0Ymgi0tYWwaSkPo9StPkYMV8NJlcheqgRIQS2/KmoZuNWCjMyMro4HrajqirZ2dk0NMQrwyiKkvII+c+cAEHvEdQIQXTHm5gOfYLQNJRYmAlmG/rIKcSEDoqKe8bnWNFLXuhURFAvmllERa2P5ngFZGxWE+cuGtn3SUMMb+QI79Tci04ME1Y0oqiY0Im1HaEAghzreE7Lvw2LmpqQmmhzOf5Dr7dtdeRIjrVWkjHxwpTY0BMul4s5c+ZQV1eXmAPy+XwUFBSQlZWFzWajpaWF/Px8KUBJRVHQZ56DPvMcoy3plQyXlasumcaGDfEl5csum2OwRR34oiEq/F6CWpR7Nr/A1yYtZqQ79elrP2n8Y0JsNOL+UR3iA7QFzjZGDnDM/zbjMpamxC4t0PPQX/PXxRc4DMgF3Y7H4+k182FhYSGFhcbUrTd+bVfSBV8gQjAURdcFui4IRWJp44T49KEt+GMR9DZ/oF/tessQO05kaJXKYZg5swSU9hXUDrGxeEYaKj7pzNDqAaU5G7dVsvGTShQFTP74Ktivn9qGrgtmTszn/MXd82mnkpZI12Tmx/sEpYr5uf+Ht2ruIqy3Ylc9bSMuQUhvBgQmrJhVB/n2aYx0JW/B4HjMrnwyJl1MLFCHyVmAHmoCBBbP6JTZMNgYcgKkeGtQ6ssRVgdKQwVKcx36qOmIMcbUxm5H03Q2bYsHyHaeK2yfgN6+r46FM4alNCHZ8SwbMZXnP96XyAz0hVHGpHBwmHNYNvxhQ+7dHyZHNiZHW50th8dQW3qi3R9I0zQyMzPJzMykrq4Ov9+fyA+dyt7akBIgpeYIln8+hKJFEXR0kkXpRmKnXYI++3zDbFNVBZvVTCgc63G/yaQYXpxwsmcYOzPziegxrl94cdpEwqcjMV81WqgZS9YIVEv8OWkhLzFfDWZXPiZH6ufOgsEgH374YaIeGMRXyFpbWxPbBQUFTJ+eOpePISVA6uFtKFrcH6KzxiuAae8HhgqQoihccu443vm4HF1AS5UJBBQWutE0wWkzi7DbjH+7VEXBbrJI8emDSONBAkffBiBU7SRj8kpENEDrvn/Fk5EpKu4JKzC7UptipaGhoYv4AF3EB6C2tjalE+bGf6JTiOgjcbrIN36pe3hBBqtXTAFgw4Y9AFy2zLgYK8nJEW0+lvhfRANogXq0kDcuPgBCJ9ZSkXIBysjoHjJjsVi6OCm63W45BEsW+vi5RHUNteogKAo0lKMEWhFF44mdc4XR5vVKbWMAb2uYUUUZ2KypfctqAi3sa46n7pzoKSAmdLzhAJtqDnNa4ZiU2tKZI763qA5uxW0aRnP0KDY1C7PJhoqZsO5HBaZ51mA3pz6LgMlVQNR7JL6hWjA5slFMFtr9k+LHpD7BXFZWFjNnzqSqqopYLEZubi4FBQUcPXqU1tZWMjMzeywcmkyGlAAB6BMXoE9cYLQZA2bfkUZeePsQQkBOlp0rLpySsrmggy11/HL763EnTeI+G6NbWhDAH/dt5EBLLVdOGFiZ7lPJ9sbHOOx/DYDuSUw6qAh+yNLiB7GZUhsVby+YhmqyooVbsGaPiYdjWJy4xl9ArLUSs2sYlkxj4ufy8vK6FSKcNGmSIbbAEBSgwUbp4cbEqlhjc4iahgAjhqUm+nxLfVlCfAB06FIbY0t9eVIFqLccxw3hRgSzBnSNF80v9ihAycpx3I41d0K3NktGsWGVUdMV6YiY5uTndEz2WswqnozULcOPaMv73BsFKUpbezzmE3JETK3bgh4LEW2tRI8G0bUYwZpdBKu3oQWbUmrHYGFo9YCEjnJ0F6Y9GxFaFL1gNKb6MnTPMMjKQ+SPQBQY6+x3PItmFGExq3hbwkwdn0uGK3WevacVjiEmdDbXHQNFMD9vFPsrP6Ap5GdadhHXT0muk19vOY41Pcp7dT+lNVqOqliI6SEUVEyKFUVR0UQYVTEzNevLjHKf1eO1k5HjWI/4aN37L0QsiGKyIYQOenyCN1y1BceIJdjyjBvudKapqQld18nJyTHUS3tICZD5jccx7f+oo6G8FADTsV2JpugZq9Cn9/yhNQJVVZg/LXXpEY7njGHjOGPYuMR2jWUbmRY7l003Lp7OpFo4q/BOw+7fG1HvUUQsCBxXB6yNSH1pWgjQgQMHOHYsvlJXWFjItGnG5KmGoTQE03XUzuLTC6Z9/R+TLILhGK9tPMIjT2yloqaVylofL759UKbjGCSotr5X3PrbnyqqqqoS/9fU1PSafiYVDJ0ekKoisoehNPVdKUHk9Jey7OToK1FaO5X+QiJ6fM7CTDwJfOnhJo4ee4dCZ32f5yY7UZqkfyxZJThGnE6stQKzuxBFtRKo/Ah0DYtnFI7hC402EYj7+jQ1xeekXC6XYWWZYSgJEBC98AbM7zyDWr4HhEBYbCixKMJig6w89OGT0OYmJ9NfX4nSIB7/1S4+ALqlozBiWLP1+SuVikRpzx3+hNcr9zKiOf7Bvendp0CJ14U/u2gCl4+bl3QbBgO2vInY8iYmtq254w20pmemT5/OkSNH0HU95X4/xzOkBAh3NrHl/23Y7XtNlNbGv94+SnmNHwDd0uGoNmFUDucs7D3w89MkShtIz8xr1vkgO54SpH0wGBVaYuONir1Etx/GKbr/kiarZ1Yd3Epl4GM81tGMzfj8Kb/+qURoUYQeTcSEGY3FYmHChO5uAkaQdAF6+umn2bBhAxaLhfvvv58RI+LhEC0tLdx8881Eo1GEENx5552GToalA8tPH8HuQ02UHvYSjWpYrSbGlWQxZ0pe/yefJP31zACinYr+NeX0sKytxI85vpOWrJ6ZN3KED+rXA4KywLuoioXR7s8l5V6flmhrFf5Dr4Eew5ozAeeoM4w2Ka1IqgB5vV6eeeYZnnjiCXbv3s26detYv349EM/I/7Of/YzCwkIOHjzIj3/8Y/74xz8m05y0x2RSmTEhlxkTclN63/56ZroQtO5/j/2t9Xh7EKCxrlwun9e9l5OsFLat0Uo6u0S2RsuTcp9TQbh6WyIGLNK4H1vhDEz29JiMTgeSKkDbt29n4cKFmM1mZs6cyeHDhxP77HZ7ovyH1WrFZEpdqgmlYl9bTiA76rE9YLOjT1yIKDa2W7ptXwMf7qxF0+JfLgWwWVWyMmxkZ9hYPKsQmzX1KTlUReGGiWfgj0VQgJjQsatmBBDWYmRaU1tqJt8+Dbsph5DWiEmxMtyZ+nCQgaJ0TkavqG0xYcbg9/upra0lOzubYDCIz+cjOzsbVVWJxWJYLBY8Hs9nJxi1ubmZrKwOtT8+K3972wMPPMC11157QtfeuXNnj+0+nw+ns/extum9ZzHveDN+bzql5SjdRGzO59EW9V3O1+fzsXnz5hOydSB21XtDbNzWNbJJAKGITqghSE1DEE0XnLeo5xiiZNnVGZe5uxOkzdT3R+hk7erLNrspi88V/oimyAEyLMW4zCeXzzgVz8xRsgiha4ioH1vhzAHNAyXDrqamJrZu3QrQpSNQVlbW5bj8/HxmzOh5vvHTvJe9kVQByszMZO/ejnrdPS333XfffSxcuJDTTjvthK49ffp0bLbuw4E9e/bQ0NDQa/dfqapEMfe8MiF2HUM09D5saI8fmjfvxFd8jh492udKli/QvW7T8bQGes8N7Xa7k2LXp+Vk7YK+bbOZMhjm+HQJ+1PxzFSLE/e4E8szlQy7yssHNkytq6sjEolgtXb/senNrnA43GuHoD+SKkCzZs3ikUceQdM0SktLuy35/frXv8ZkMnH11Vcn04yumC3QWy7jHn7hU8XwAhcOu5lgqOeMiKpCyueGID7E+p89b9IUDXJG3hjOHjaeFyp20xQJ4rHa0YXgjPwxjHan3jbJwMnIyOh3tRPi0yFmc+oWx5N6J4/Hw8qVK1m7di1ms5n77ruPDRs2UFJSwogRI1i/fj3z5s3jK1/5CgUFBfziF7/41PfsLX4ogb8Z8ztPoTRUQiSMEvaDakIvHk/svDXQR53zZMQPtWMxq1y5YjwHy1uobQhgt5lw2ExkZdixWU04bGbcztTPH9y57cVERPwbtQd4s/YAx//Gbm4s55uTz2ZUP8GrEuMYPXo0wWCQxsZG7HY7sViMcDiMw+HAZrMRDodxOp2MHj06pY6JSZe6NWvWsGbNmsR2517Qnj17kn377riyiC0zzheoL0wmlYmjPEwc5THaFABCWrRLOg6gm/i0s9NbKQUozZkyZYrRJnRjaDkipjHb9zXwftsktNWiYjGp+EMx7FaVghwntY1BSoa5OGfBcExqalYp7CZLpxx+fTM+I7UZ/oQQbG/6M5XBj/BYxxCMNdEa65jnMCt2Tsv7Nrn2iX1cRWI0UoBSRG/JtdqprPMn3gwdCBN/c2IBqPTG24944bkKG05717ctmcm1fjD9An6+5w1CWowJ7jwuLpnOhvLttIRDuC02VEXhjIKxTMos6P9ip5Dq4BaO+P8DQG1oR7f9MRHi48ZHWFr8YErtkpwYQ0uAdA2l5ijCmYESCYK/GWwuROEoUI0teTNQenJlSCbZNif3z76oS9s3JhmfrkRH6/cYIfo/ZqghhKClpQWz2YzD4aC5uZloNIrb7R6wa8GpZOgIkK5jefHXqOWlCBSUTgMLfcQUosuvhyROvvU3Of7Bzhq27mkA4r5JFrNKJKZjUuO5oOuaQgwvcLH8jBGYTV3tTObkeLpS5JhHkWMeVcEtZFlGEdKaCesdWQcVTMzOvsZAC9OT3bt3U1MTH+rb7XZCoY5qt1OmTKGoKDnZIHpjyAiQ4q1GbUtAphw3q6GW7UFprkVkG5f4a9H0QhZOK0DXdVRVRVEUYrqOuU1YdCFQZX3xBKpiYmHeNxBCR1Hiz0jTIqiqBSE0VNXYj3aoZgeRxoMA6CEvIFBtWbjGfR6TLXk5vRsbG/sMgWmPPvD5fPh8PtzujrS6u3btSjgr9nTd45PZnwqGjAAJR2Y8/UYPPkDCYkc4U1s54XhC4Ri6iE9AN/siKIqC22khGtPx+SOEIhr5OY5uvZ9kUh/288i+d2mKBLvtUwCrasZuMrNm9NyUzwEl7FA6nofJZG1rS/7Huq8veoYlwoy8eG/saEP8x25UroIebqZ86/Psauw98PfTfNHz8/tfCGgvOthekLCzAAkheu1J5+XlDej6J8qQESAcbqIrrse0/c24EEVCKN5qREYu2txlYGClz9LDTby1uQohwGxSiLXFgmW4LASCUbS2tW+zSWHNsnG4nKfOYbKvyfHqYAueWARPP9d4Y++/Oejq/qVKduUJo+jvi9i5atLRtoQAo9r8NC1q719y+HRf9IGkPWlpaWH//v00NTVhsVjIyckhFAqRm5vL1KlTsVhS62s2dAQI4gUIi9IvQdSHO+sSpXfaxQeg1d81PCOmCT7cWcc5C42pKdUbKZ4XN5z+vuhC12jc8wJhb9lxexRKZi5lfM7opNnWH5mZmcybN4+jR48CsGTJEsNsgSEmQOmK3WYi0EsIxvE4Haf2Letrcrw25ON/975DS6zn0BWLomJVzawdM48pWd0DQpM5OR7RfNSF92BTM7CpWTRHywjrXpymfKyqi4BWT4F9BjZTamqodUZRTeRMvRg9EsBa8zJC18ibeTZmZzaqybhwn3REClAacP6i4by7tZqYpmO3mqhtCqEqCmNLMmjxhamoDSKEoDDXwcLpqZtrKbC7uWfW8pTdb6C0Rqt4q+ZuNBHq8zhVsXB2wd1kWktSZFkHiqJgsrkABUU1Y804uYj9zzpDS4CEQN3xNgKBGosgYhFEbjGEQ+DKArsLkV0IVkdKzcrJsvOFz41O6T0HCz1N9mqFO9ALuopPuDL+BbcVd6Q00UWUt3b+BVP17B6vm4xVnVDjEQI1e0BRELqOFmrBZDd2gQM6/H+CwSDRaNehva7rtLS0YLPZcDhS+9kfUgJk/fP/RQn5+jxGuDxELr0V3KmNawqEYlTU+qiuDwKC4nwXqqqgaTq+YIysDCsjh2WkLAwjHehtMlZE3d1i0kKVcReKzgIEoMbcPQ4Dk7GqEwt6aSx9CTrFz2mR9Jgg27p1K5WVlYntpqYmsrOz0XWdDz/8kPr6elRVZf78+RQUpK6X/ZkUoJ5+NTMiLZzfg/jsVeOrN5P0+HKF4vey/19/ZX/m2B6vm4xfTV8gyt9fPUgo0vHB3XXQ2+24EcNcrDhjZMoy1oW1GL87sInKYDNFjkwyzDYaIwGiuk6Gxcqy4imMSWIajt4me4UQbK9/giMtb2FWndhNHloxIdCxmTIxYUMnysiMJcydeFWXpfpkooVbu4hPwl59YPN7ySIYDHYRH4j3esLhMC0tLdTX1yfajhw5IgXo09Dbr1rE4uiaAbGN0uMECCBg6blWUrJ8IcqqfV3Ep/fj/ITCGg57at62n+56naZo3AfooK+hy76qEOzf+w7/d/rnybOldqldURRm5V/BrPwrEm3Njg0AXDb++ym1pTPWjCIsrnyi/q55d1Srsa4IFosFi8XSbehlMpmw2WyoqppIZCaHYJ+SvpZItS358M7fQW+LERIC2n8dbS7IGYYyYS6L5y9LgaUd5GT1UGmiB5x2E9YU5oT2Rrs7IHZGAMd8jSkXoHRFMZnJm3EZEV8tJqsDLRLEUvUeisFxhmazmYULF7J37168Xi+KouB0OhPxYPPmzePIkSM4nc6Up+z4zAlQX5jmng9zj0uPuSH+y2m67A4DLIpTmOtk6ZIR7DrQSENzfHI1N8uGzWomFNEIR2J4MmzMn5af0jmgXJuL+rC/1/0qisyEeByKyYwtqxgAsyPbcPFpJzs7u0va4w1tn3uI14cvLDRmlW5ICVA6M2Z4BmOGp95npS/umHoejx36iKOBJka4ssg0O6gL+4jqGpkWOxcUTSbHQA/ydEToGv6qHURaq8FsIxZoROgaTXtfxWRz4cifiMWVvDpvgw0pQGnE4YoWPtpZSyCkIRAU5Tk5d+FwKusC1DYGGVnkZlhu6r7wZlXla+PTs+RNMNrEhzWPEtWCeOyj8EXjJYOOtrzHqMzTDbOr+dDbBGp2J7b1WHwVLFi/DwB/9S4K5l6ByeB5oXRhyAqQiAQR3rpEgnr9wFaU4RNQHO5+zjx5+gpgDGsWakLFtE+Tq9E6jvnhT5VN6MTjc7bsqWOYoxKr2nUyMVmrc+mKEIIXj3ybiB4PqKwL7SYQmwXA+1UPEtWDjPecWCWKU0WktabP/UKLEAs0GSpAgUCAUCiEyWQiHI5//iORCG63O6U1wWCICpBorkN/4gEINENbiR7xz9cQior65dtRhp/6AoX9rZ6Fo046r9Gp0VoAYpbOwqIQ1e3YzV0TbSVrda4dXQhCWgREPA2GzWRCURQ0oaOgpDxNSFQPJMSnJ+qCpYYJkD13DL5AQ6/7TTY3Fndq09d2pqGhgY0bN+L3x+f2Xn311cQqWH5+PgsXLvzsFCZMV8TO9+Li022Hjv7BvzBd9q1Tfs/+AhgbvAH+/Pzubu3ZmXa8rWFEW6qOL6w4m+zM1FUi3VR/hKeOftKlzYTC5wrH82btAcyKiSvHzGO6J3WJrKwmF05zHoFYfQ97FYa75qbMluPJHLkIq7uQqL8exWTDdGwLQtfIHDsfVTVjyx6Nah7Yqmcy6Kl+V/sSfF1dHc3NzXg8npTZMyQFiAxP7/uyjBnK5HqcXLFiCpu2VdLsCxMIm7DbTFx1yXRqGvzUNAQYVZSZFPHpa2j4Rn4EOrlEeRrjXfbXxX5QQBMx/rrvA86s7x5kmcyh4YVjHuST2seJiSAe22g2W+JDn/NGXEaBc2pS7jlQ7DmjsbdFvJts+wFwF8000KIOnE5nIhfQ8SiK0mOxz2QyJAVImXEWtDQgDn4CPkvcH8jqgOHjUc9cZZhdw/JdrDw/PvzbsCE+aamqCkX5borykzM31d/QTYUu2Zez2wTIm93xQTUJJaWOmwBm1cb8YR0pV3eb48vKRotPNNBEoGYXismKs3AaWrgVoUVoKduMq2BiUrMhDoT58+fzn//Ek/krikJhYWFiLmj06NHSETEVKIqKcsYX4YwvdvID+o7BVsXRdcE7m8uoqvOjqvDEC3uoawpiNissmlHEvGmnNm1sf0PDSY2V/O+ut9A6pbE1qypXTljIy+W7sCgmvjrxNMZmDp1J8N6IBhqo++TpRDiGr3wzWjgu375jm/BXfkLBnNWGTkArisK5556L1+sFYMGCBYbZAkNUgNKZLbtr2Ly7FrOuo+lQVR+fLIxp8NbH5eR6HIwenpUye6blFPPImR2FJTfUxQX7zKLxnJmGyd2MJOwt7xoLdlxcmIiFiLRU48gbl2LL0pehVUphENDs66VufRstvkiKLJGcKFZ3P0GcigmLS3qOd2ZI9oD0hgpEXXl87ifQEq8NX1+OYnehpDgNx/FMn5DHrgM9L+M6bCbGj/Sk1qA05q3yB6gJ7MBhziGsLSAQq+f1Y3czv/BasmypT0JmzSwiZ+rF+Cu3oZgsuIfPQz38PEKPYs8bh7t4FmaHJ+V2tRMOh4lGo3EXCk1DCEE0Gk15HujODDkB0v79J9j5TkdDux/QX+6KR8ufuQp1QWqDUTtTmOvi2i/N4LnnDqCqCsuXT+dwRTNuh5VRwzOxWlIbW/RBzWGeOrQZTeiMdGZj9XsB+O2ed5mdW8KCgtEptaed/U2vUOnfAoAvWkNzJJ5/uTa4i41VD7Fs9M8MscuePRJ79sjEttkZz7aQM+kCQ+xpp7y8nG3btiUKW7bPAb366qssWrSI3FxjemZDToDY9V6fu8X7/wADBQjAabfgsMXfGk+mnTkp9PvpTLm/iT/s25jY3tdax5hYfAi4s/4YH9cfw6SqzM0b2dslkkYg1ruzX1jr3Ukx2QgtRsRXg8mWgaKoiFgYxWQl6q9HaDGsmcbUntu3b1+PVXV1XefgwYNSgFKG2ZIIv+gRV+rTZ9Y1BXjpncOEwjFOm1XMxk8qCTfEv0S//fs2rGYToYjG6XOGM31C6labDjTX9XvM0dZGQwRoSs5K9nv/TVT3Awp2UxYhrRkFEzPzVqfcHoiLT/2ODW35gFRQFaKBGCgqdZ88BYCzYAqeCeem3DabzUYgEOh1n1EMuUlo5eIbwdJDj8Jig2FjUC++MeU2vb7pGPVNQXyBKK9tPIo/2BHr1eqP0tAcwh+M8urGIwQHWD3jVDA7t6RbArfOmBWVmbnGlAiympxcOu53fH7k/Vw27o9kWoeTZ5/EZeN/x5issw2xKdJa3SkZmd4p71THaligdg9CS32GxNmzZ5Ofn4/NZsNut6OqKqqqUlhYmPIcQJ0Zcj0gdfQ0uPl/OxoSfkDfM8gi0LT+syFCfM5cT2ERLo/NyY/mX8yTBz4iomvMzBnOsfItKMDicQuY6CmgyJk6l4DjMalm8hwdcXuqYsJqSl4wcb/22DPiCe56SMuaOMbqBgNyBLlcLhYt6shs0J4PSPoBSTh7/gj++eZBIhGNhTML2byrNuH2ZzGrqKpCJKqxZPZwXI7UrlgUODL4xoyOIcMGW2nc5uJTH7A72DHbs8iZsoJAzR5M9iwU1YR6eDuq2Y6zcDhCi+IesSDlEefpjBSgNKBkWAb/5/JZCBEPvVgyu4Rnnz0MKHzxi3MRQiT2SdIbe/Yo7NmjEttmZzkAnvHnGGVSWpP0OaCnn36a1atX85WvfIWysq6lardv387q1au5/PLLE/EpyUZEI4jWprYvtQAtBrqG7q1Fb2nscaUgFSiKgqoqNLUGKT1c39bWdZ9R+KNhvOEATeEAmtAJazFKm6qp8HsNe14QX+0Ka634o3XERFteJ6HRHK7AH+1/Aj1VCC2KFg0hhEAL+xC61v9JSSQWi6HrOpqmUV9fn6iQcXzS+lSQ1B6Q1+vlmWee4YknnmD37t2sW7eO9evXJ/Y/8MADrF+/Hrfbzdq1aznrrLMwmZI3PhZ1Zeh//wUEW2HsrPikSlV8ZUD84eX435FTUS+9BcWU+s7hy+8cZveh+PKyOeCjIMf4dKf/qdzLkwc3J7bHtLQA8ObONwAodmZx59zlmFJU+qad0sYX2Fr3Z2gbrDaHZmFRHLx4+Fu0RqvitmZ+jtOKUr+o0I4W9hNpqQIENR/+HpPDgxb0YrK6yZ1xKWYDChbW1dXx8ccf09TUBMCmTZtQFAUhBIqiMHv2bIYPT93CQlK/Zdu3b2fhwoWYzWZmzpzJ4cOHE/vC4TCapiWSYY8ePZojR44wblxy4mRKS0vZ/e7rEC0EcyEci8dY1Svx6N/n2xwSqYzAU08wde4CJk+enBRbutm1ezdCQEVta+INUbQA9fWBLsnDAaZOnZpSuw601DGmUy/HEYyv4Iw50NLW0sITB5/GabamxLZ2u+qCexB0pLiItbqJAeUbC4B4SMQnNFFuf4rp02alzK7OaKEWvIH4s3t7vwCa2va0Ytr/LDPmLk7pewnQ0tJCNBolEon7c1VVVXU5tqamhiVLlqTELkiyADU3N5OV1bFK0rm77vV6ycjoSE2QmZlJc3MPScJ6oafESn1RV1dHRNNpz1rTXiPMSfcl0UAozJEjRxJZ45JJXV0dPp+v+1BGsaCo4PN1LaaYaruU48yKWrr3dMLBELoSSYlt7XahqqB0DGVUaw8xcgL8/kBq7eqEqkex97JmEInGUv5eAmha/Jn1NtIQQqTMLkiyAGVmZrJ3797EduecMVlZWV0SI7W2tnYRq/6YPn36CTtQifBZiFf/gmiqjucEAsS2NyHsjzsnmm0o85eizjPGbf7gsSb++eZBhAC7zcR1X5yJJcWhF8dzpLWBJw58RH3Ij0U1EdVj+GMRBPHMiOeXTOayMXNSbld9cD9bav9IVA8QiDUihM60nMsIak0cbnkTFQsLCq9jZObilNvWjtA1GktfJOytwGR1Ys8dS9hbjsWdj2fc2Shq6of5wWCQ7du34/P5CAbjdd8yMzNRFAWTycSMGTO6dAwGQjgcPuEOQTtJfQKzZs3ikUceQdM0SktLGTWqY3XAbrdjMpmora3F7XZz9OjRLvuTgWJzolx0fdfG2an3Su2NcSOz+dZX5xttRhdGZ+TyvTnGhqb0RJ5jAheMur/HffMLr+mxPdUoqoncqRcbbUYXHA5HF38go0mqAHk8HlauXMnatWsxm83cd999bNiwgZKSEhYuXMh3v/tdvvGNbyCE4MYbb8Rsll4BEslQQhFGrqOeBO3dvZMZgkkkklPPp/lODrlYMIlEkj5IAZJIJIYhBUgikRiGFCCJRGIYg27ZqX3OvN2TUyKRGEv7d/Fk1rMGnQC1B8zt27fPYEskEklnotEodvuJpQ8edMvwuq7j9/uxWCwyr4pEkga0V9dwuVw9Vsjti0EnQBKJ5LODnISWSCSGIQVIIpEYhhQgiURiGFKAJBKJYUgBkkgkhiEFSCKRGIYUIIlEYhhSgCQSiWFIAZJIJIYhBUgikRiGFCCJRGIYUoAkvVJTU8Ptt9/e6/7y8nJefvnlpNqwZ88e3n///cT2+vXr2bJlS1LvKUkdMhhVctJ88MEHPPnkk/zP//zPgM/RNO2Eym9v2LCBQ4cOcdttt52MiZI0RwqQpFfKy8u59dZbWb16Ne+88w6NjY1UVFRw4403cumll7J69WoOHjxIcXEx11xzDeeddx533303hw4dAuAHP/gBs2fP5rvf/S52u50dO3Zw4YUXMnLkSB599FEikQjFxcWsW7cOl8tFTU0NP/jBD6iursZsNvPwww+zdu1aIpEI+fn5fO973+Mf//gHK1as4KyzzuKtt95i3bp1CCFYunQpN998MwCnn346y5cvZ+PGjYwcOZJf/epXJyR6khQiJJJeKCsrE6tWrRLPPvusuOiii0QgEBC1tbXinHPOEUIIsWnTJvHNb34zcfzPf/5z8eqrrwohhKiqqhKXXnqpEEKIO+64Q9x6661C13UhhBBerzdxzm9+8xvxpz/9SQghxM033yyeffZZIYQQwWBQBINB8eyzz4qf//zniePvuOMO8dZbb4lgMCjOOeccUVlZKSKRiLj88svFRx99JIQQYuLEieLDDz8UQgjx9a9/Xbz77rtJeT6ST8+gy4goMYbFixfjcDhwOBzoup7ITNmZ999/n7fffpuHH34YAK/XSywWA2Dp0qWJBHJVVVXccsstNDQ0EAwGWbJkCQBbt25l/fr1AP1m1jt8+DDjxo2jqKgIgBUrVrBlyxbmz59PZmYmCxYsAGDKlClUVFScgicgSQZSgCQDwmq1Jv5XVRVd17sdI4Tgt7/9LYWFhd32dRaUH//4x9x8880sWrSIl19+mTfffDOptmqadkqvLzl1yFUwyUnjcrnw+/2J7cWLF/O3v/0tsV1aWtrjeT6fj/z8fHRd5/nnn0+0z507l+eeew6IV9sMhULd7tHOmDFjOHjwIDU1NcRiMV5++WXmzp17ql6aJEVIAZKcNJMmTSIUCnHJJZfwz3/+kxtvvJHa2louvvhiVqxYwTPPPNPjeTfccAPXXXcdq1atoqSkJNH+/e9/n5deeomLL76YNWvW4PV6WbRoEdu3b2flypVs2rQpcazdbueuu+7iuuuuY+XKlSxevJj58+cn/TVLTi1yFUwikRiG7AFJJBLDkAIkkUgMQwqQRCIxDClAEonEMKQASSQSw5ACJJFIDEMKkEQiMQwpQBKJxDD+f9rfIFH2MvH8AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "# Plot the boxplot\n",
    "\n",
    "sns.set_theme(\n",
    "            style='whitegrid',\n",
    "            font_scale=0.8\n",
    "            )\n",
    "sns.set_palette(palette=sns.color_palette([region_colors_dict[r] for r in interact_with]))\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(4,4))\n",
    "sns.stripplot(result_df,x='interaction',y='active_in',size=4,hue='interaction',ax=ax,legend=False)\n",
    "sns.boxplot(showmeans=False,\n",
    "        meanline=False,\n",
    "        meanprops={'color': 'k', 'ls': '-', 'lw': 2},\n",
    "        medianprops={'visible': True},\n",
    "        whiskerprops={'visible': True},\n",
    "        zorder=10,\n",
    "        x=\"interaction\",\n",
    "        y=\"active_in\",\n",
    "        data=result_df,\n",
    "        showfliers=False,\n",
    "        showbox=True,\n",
    "        showcaps=True,\n",
    "        color='whitesmoke',\n",
    "        width=0.6,\n",
    "        ax=ax)\n",
    "\n",
    "ax.set_ylim(-0.05,1.05)\n",
    "ax.set_xticklabels('')\n",
    "ax.set_title(title)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig('./plots/receptor_ligand_interaction_analysis/Club_epithelium_as_target_ligrec_pct.pdf')\n",
    "plt.show()\n",
    "\n",
    "result_df['direction'] = 'club_as_target'\n",
    "plot_df = pd.concat([plot_df,result_df],axis=0)\n",
    "plot_df.to_excel('./source_data/figure_4g.xlsx')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "result_df['n_interfaces_with_club'] = result_df['interaction'].map(n_interfaces)  # Add a column with info on the number of interfaces\n",
    "\n",
    "# Calculate p-values for overrepresentation of individual ligand-receptor interactions\n",
    "list_of_pvals = []\n",
    "for i, row in result_df.iterrows():\n",
    "    df = result_df[(result_df['ligand'] == row['ligand']) & (result_df['receptor'] == row['receptor'])]\n",
    "    df = df[~(df['interaction'] ==row['interaction'])]\n",
    "\n",
    "    a = row['active_in'] * row['n_interfaces_with_club']\n",
    "    b = row['n_interfaces_with_club'] - a\n",
    "    c = (df['active_in'] * df['n_interfaces_with_club']).sum()\n",
    "    d = df['n_interfaces_with_club'].sum() - c\n",
    "\n",
    "    arr = np.array([[a,b],[c,d]])\n",
    "    arr\n",
    "\n",
    "    stat, pval = fisher_exact(arr,alternative='greater')\n",
    "\n",
    "    list_of_pvals.append(pval)\n",
    "\n",
    "result_df['overrep_test_adj_pval'] = multipletests(list_of_pvals,method='fdr_bh')[1]\n",
    "\n",
    "result_df = result_df.sort_values('overrep_test_adj_pval',ascending=True).reset_index(drop=True)\n",
    "\n",
    "dict_with_final_ligrec_results['club_as_target'] = result_df\n",
    "\n",
    "# Save the results\n",
    "result_df.to_excel('./source_data/club_as_target_ligand_receptor_activity.xlsx')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Kruskal-Wallis H-test: test statistic = 59.2094, p-value = 0.0000\n",
      "             Category1           Category2      stat       p-value\n",
      "0           Fibroblast    Basal epithelium  5.111546  3.195330e-07\n",
      "1           Fibroblast  Luminal epithelium  5.027397  4.971817e-07\n",
      "2               Muscle    Basal epithelium  4.923679  8.493209e-07\n",
      "3   Luminal epithelium              Muscle -4.344423  1.396424e-05\n",
      "4           Fibroblast               Tumor  4.272016  1.937140e-05\n",
      "5               Immune    Basal epithelium  4.203523  2.627930e-05\n",
      "6   Luminal epithelium              Immune -3.808219  1.399712e-04\n",
      "7               Muscle               Tumor  3.616438  2.986844e-04\n",
      "8               Immune               Tumor  2.788650  5.292828e-03\n",
      "9     Basal epithelium               Tumor -2.146771  3.181151e-02\n",
      "10          Fibroblast              Immune  1.780822  7.494155e-02\n",
      "11  Luminal epithelium               Tumor -1.414873  1.571058e-01\n",
      "12              Immune              Muscle -1.107632  2.680207e-01\n",
      "13          Fibroblast              Muscle  0.743640  4.570943e-01\n",
      "14  Luminal epithelium    Basal epithelium  0.655577  5.120961e-01\n"
     ]
    }
   ],
   "source": [
    "# Perform Kruskal-Wallis H-test\n",
    "interaction_categories = result_df['interaction'].unique()\n",
    "grouped_data = [result_df[result_df['interaction'] == category]['active_in'] for category in interaction_categories]\n",
    "test_statistic, p_value = kruskal(*grouped_data)\n",
    "print(f\"Kruskal-Wallis H-test: test statistic = {test_statistic:.4f}, p-value = {p_value:.4f}\")\n",
    "\n",
    "\n",
    "# Perform pairwise Wilcoxon rank-sums tests\n",
    "interaction_categories = result_df['interaction'].unique()\n",
    "pairwise_tests = list(combinations(interaction_categories, 2))\n",
    "\n",
    "# Collect p-values and test statistics\n",
    "test_results = []\n",
    "for category1, category2 in pairwise_tests:\n",
    "    group1 = result_df[result_df['interaction'] == category1]['active_in']\n",
    "    group2 = result_df[result_df['interaction'] == category2]['active_in']\n",
    "    test_statistic, p_value = ranksums(group1, group2)\n",
    "    test_results.append({'Category1': category1, 'Category2': category2, 'stat': test_statistic, 'p-value': p_value})\n",
    "\n",
    "# Create a dataframe from the test results\n",
    "test_results = pd.DataFrame(test_results).sort_values('p-value').reset_index(drop=True)\n",
    "print(test_results)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save the results\n",
    "pd.concat(club_as_target_ligrec).to_csv('./plots/receptor_ligand_interaction_analysis/{}_as_target_ligand_receptor_results.csv'.format(source.replace(' ','_')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Save both directions into a single excel table (Supplementary Table S9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a Pandas Excel writer using XlsxWriter as the engine\n",
    "writer = pd.ExcelWriter('./supplementary_tables/club_regions_ligand_receptor_analysis.xlsx', engine='xlsxwriter')\n",
    "\n",
    "# Write each dataframe to a different worksheet\n",
    "for sheet_name, df in dict_with_final_ligrec_results.items():\n",
    "    df.to_excel(writer, sheet_name=sheet_name, index=False)\n",
    "\n",
    "# Close the Pandas Excel writer and output the Excel file\n",
    "writer.save()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chord diagrams (Supplementary Figure S8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#from d3blocks import D3Blocks\n",
    "\n",
    "# Get Club region degs\n",
    "club_markers = check_top_markers('Club epithelium')['gene'].tolist()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Club as source "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download from Supplementary Table S9\n",
    "club_as_source = pd.read_excel('./supplementary_tables/club_regions_ligand_receptor_analysis.xlsx',sheet_name='club_as_source') \n",
    "\n",
    "club_as_source = club_as_source[club_as_source['ligand'].str.contains('CXCL|CCL')]\n",
    "club_as_source = club_as_source.sort_values(['ligand','receptor'])\n",
    "\n",
    "# Add the weight parameter ()\n",
    "club_as_source['weight'] = club_as_source['active_in'] * club_as_source['n_interfaces_with_club']\n",
    "\n",
    "# Concatenate into a single dataframe for plotting into a single Chord diagram\n",
    "plot_df = club_as_source.groupby(['ligand','receptor']).sum().reset_index(['ligand','receptor'])[['ligand','receptor','weight']]\n",
    "plot_df = plot_df.rename(columns={'ligand':'source','receptor':'target'})\n",
    "\n",
    "plot_df.to_excel('./source_data/figure_4e.xlsx')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot\n",
    "# Needs access to a browser due to how D3Blocks is built\n",
    "\n",
    "d3 = D3Blocks()\n",
    "d3.chord(plot_df,save_button=True,\n",
    "color='source',\n",
    "ordering=plot_df['source'].unique().tolist() + plot_df['target'].unique().tolist(),\n",
    "title='From club to all',\n",
    "fontsize=30,\n",
    "cmap='viridis',\n",
    "arrowhead=50\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Club as target"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download from Supplementary Table S9\n",
    "club_as_target = pd.read_excel('./supplementary_tables/club_regions_ligand_receptor_analysis.xlsx',sheet_name='club_as_target')\n",
    "\n",
    "# Plot the receptors that are within the top 100 of Club region markers\n",
    "club_as_target = club_as_target[club_as_target['receptor'].isin(club_markers[:100])]\n",
    "club_as_target = club_as_target.sort_values(['ligand','receptor'])\n",
    "\n",
    "# Add the weight parameter ()\n",
    "club_as_target['weight'] = club_as_target['active_in'] * club_as_target['n_interfaces_with_club']\n",
    "\n",
    "# Concatenate into a single dataframe for plotting into a single Chord diagram\n",
    "plot_df = club_as_target.groupby(['ligand','receptor']).sum().reset_index(['ligand','receptor'])[['ligand','receptor','weight']]\n",
    "plot_df = plot_df.rename(columns={'ligand':'source','receptor':'target'})\n",
    "\n",
    "plot_df.to_excel('./source_data/figure_4f.xlsx')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot\n",
    "# Needs access to a browser due to how D3Blocks is built\n",
    "\n",
    "d3 = D3Blocks()\n",
    "d3.chord(plot_df,save_button=True,\n",
    "color='target',\n",
    "ordering=plot_df['source'].unique().tolist() + plot_df['target'].unique().tolist(),\n",
    "title='From all to club',\n",
    "fontsize=30,\n",
    "cmap='viridis',\n",
    "arrowhead=50\n",
    ")\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "squidpy",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.15"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
